Keywords: données transcriptomiques, estimateur à noyau, processus de renouvellement, régression inverse par tranches, tests multiples, valvométrie. v kernel density estimator, multiple testing, renewal process, sliced inverse regression, transcriptomics, valvometry. vi

Je tiens à remercier dans un premier temps mes directeurs de thèse Gilles Durrieu et Jérôme Saracco. Ce travail doit beaucoup à leur pédagogie, à leurs encouragements et à leur capacité à jalonner les projets de recherche auxquels j'ai participé. J'ai ainsi pu, grâce à eux, acquérir de nombreuses connaissances tout en conservant la capacité d'exprimer mes propres idées.

Je suis reconnaissant envers Gérard Biau et Stéphane Robin pour avoir accepté de rapporter ma thèse et pour le temps qu'ils ont passé à considérer ce document. Merci aussi à Pierrick Legrand et à Jean-Charles Massabuau pour avoir bien voulu prendre part à ce jury de thèse.

Je remercie l'ensemble des membres des équipes CQFD et ALEA d'INRIA pour l'environnement stimulant, épanouissant et productif qu'ils ont su créer. Merci aux scientifiques du laboratoire EPOC de m'avoir fourni un contexte biologique qui donne un sens pratique aux résultats de recherche présentés ici. Merci également à Stéphane Girard et à Benoît Liquet pour nos fécondes collaborations.

Enfin, plus généralement, simplement remercier ma famille, mes amis et mes collègues qui m'ont aidé durant ces trois dernières années me semblerait annoncer la fin de ce qui nous unis. Ce n'est pas mon intention.

Modélisation stochastique de grands jeux de données : applications en écologie et en génétique

Résumé

Deux parties principales composent cette thèse. La première d'entre elles est consacrée à la valvométrie, c'est-à-dire ici l'étude de la distance entre les deux parties de la coquille d'une huître au cours du temps. La valvométrie est utilisée afin de déterminer si de tels animaux sont en bonne santé, pour éventuellement tirer des conclusions sur la qualité de leur environnement. Nous considérons qu'un processus de renouvellement à quatre états sous-tend le comportement des huîtres étudiées. Afin de retrouver ce processus caché dans le signal valvométrique, nous supposons qu'une densité de probabilité reliée à ce signal est bimodale. Nous comparons donc plusieurs estimateurs qui prennent en compte ce type d'hypothèse, dont des estimateurs à noyau.

Dans un second temps, nous comparons plusieurs méthodes de régression, dans le but d'analyser des données transcriptomiques. Pour comprendre quelles variables explicatives influent sur l'expression de gènes, nous avons réalisé des tests multiples grâce au modèle linéaire FAMT. La méthode SIR peut être envisagée pour trouver des relations non-linéaires. Toutefois, elle est principalement employée lorsque la variable à expliquer est univariée. Une version multivariée de cette approche a donc été développée. Le coût d'acquisition des données transcriptomiques pouvant être élevé, la taille n des échantillons correspondants est souvent faible. C'est pourquoi, nous avons également étudié la méthode SIR lorsque n est inférieur au nombre de variables explicatives p.

Chapter 1

Introduction générale

Cette thèse contient des exemples de modélisation en biologie basés sur des outils statistiques existants ou originaux et sur des jeux de données de grande taille, principalement fournis par le laboratoire EPOC (Environnements et Paléoenvironnements Océaniques et Continentaux) de l'Université de Bordeaux. Après deux chapitres d'introduction respectivement en français et en anglais, le Chapitre 3 présente les travaux réalisés sur des mesures de l'activité d'huîtres. Le Chapitre 4 est dédié au développement et à l'utilisation de méthodes statistiques adaptées aux données transcriptomiques. Le Chapitre 5 conclut cet ouvrage. Nous y décrivons également quelques perspectives de recherche.

Un modèle stochastique pour des bioindicateurs

Dans le Chapitre 3, nous cherchons à déterminer comment résumer l'information contenue dans de grands échantillons. En effet, le signal considéré est formé de distances entre les deux parties de la coquille d'une huître, mesurées à haute fréquence. On dit alors qu'il s'agit d'un signal valvométrique. L'un des objectifs de ces enregistrements est de réaliser un diagnostic de l'état de santé de l'animal. Lorsque plusieurs animaux placés à un même endroit sont en mauvaise santé, un changement dans leur environnement est une cause privilégiée. L'étude de données valvométriques est donc une étape importante dans l'utilisation d'huîtres comme bioindicateur de la qualité de l'eau. Nous présentons plus en détail le contexte biologique relatif aux enregistrements valvométriques dans l'introduction du Chapitre 3 et dans les Comme le montre la Figure 1.1, il est possible de distinguer deux états de l'animal dans un signal valvométrique. Nous appelons ces états "ouvert" et "fermé". Lorsque l'huître est en bonne santé, les biologistes ont observé que les transitions d'un état à l'autre sont reliées à la marée, ce qui n'est pas forcément visible si l'animal présente des difficultés à s'adapter à son environnement. En utilisant un processus de renouvellement dont "ouvert" et "fermé" sont deux états, et un estimateur de sa fonction de survie, nous montrons dans A hidden renewal model for monitoring aquatic systems biosensors (voir List of works) comment il est possible de réaliser une classification des animaux à partir de leur signal valvométrique. Le résultat de cette classification s'avère être en adéquation avec l'information dont nous disposons concernant leur état de santé. Cet article constitue la Section 3.1. Le processus de renouvellement est constant par morceaux alors que le signal observé (Figure 1.1) est continu. De plus, ce signal présente des zones où il varie très rapidement vers le bas pour revenir ensuite très rapidement à une valeur élevée. Ces pics sont appelés microfermetures et leur minimum peut parfois être tellement bas qu'il est difficile de distinguer microfermetures et passages de l'état "ouvert" à l'état "fermé". Nous présentons donc dans la Section 3.1 une méthode pour estimer le processus de renouvellement caché dans le signal valvométrique observé. Cette procédure permet de déterminer un seuil. Lorsque le signal valvométrique est en dessous de ce seuil, nous décidons que le processus sous-jacent est dans l'état fermé. Sinon, il est dans l'état ouvert. Pour trouver le seuil en question, nous estimons une densité de probabilité reliée au signal valvométrique. Nous supposons que cette densité a deux modes : l'un correspond à l'état ouvert de l'huître et l'autre à l'état fermé. Nous choisissons un estimateur qui a également deux modes. La position de l'unique minimum local de cet estimateur est alors choisi comme seuil.

Plusieurs méthodes ont été conçues afin d'estimer une densité de probabilité f lorsque nous avons à notre disposition des observations indépendantes en provenance d'une même loi de densité f , ainsi que le nombre de modes de f . Nous pouvons par exemple utiliser un estimateur à noyau ou un estimateur basé sur les ensembles de niveaux de f . Dans Comparison of kernel density estimators with assumption on number of modes, nous avons étudié en simulation ces estimateurs ainsi qu'un estimateur à noyau qui n'utilise pas d'hypothèse sur le nombre de modes de f . Cet article constitue la Section 3.2.

Le nombre de microfermetures et leur amplitude peut également être informatif sur l'état de santé de l'animal. Toutefois leur détection est difficile car le signal valvométrique se révèle être bruité. Dans la Section 3.3, nous présentons une méthode existante de débruitage en utilisant les ondelettes de Haar. Celle-ci a été choisie parmi de nombreuses autres grâce aux résultats de simulation présentés dans [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF]. Nous étudions également plus amplement la détection de microfermetures directement à partir du signal bruité, dans cette section.

Pour faciliter l'utilisation de ces approches, une implémentation en R a été réalisée et une interface graphique a été crée grâce à la bibliothèque GTK+. De plus, ces programmes permettent d'importer des données valvométriques depuis un serveur SQL distant où elles sont stockées comme nous l'expliquons dans la Section 3.4. Des pistes y sont également décrites pour optimiser le temps de calcul lors d'un travail avec R.

Des méthodes statistiques pour analyser des données génétiques

Les problèmes biologiques abordés dans le Chapitre 4 concernent les relations entre des données transcriptomiques d'une part et des variables diverses que l'on nomme variables explicatives, d'autre part. Nous appelons données transcriptomiques des mesures de la quantité d'expression de marqueurs génétiques. Nous cherchons à expliquer les variations de ces expressions en utilisant les variables explicatives. Dans ce contexte de régression, la grande taille des jeux de données se caractérise soit par un grand nombre de variables à expliquer, c'est-à-dire un grand nombre de gènes, soit par un grand nombre de variables explicatives. Deux exemples de problème génétique sont fournis dans la Section 4.1.

Pour expliquer les variations dans l'expression d'un seul gène grâce aux variables explicatives, il est possible d'utiliser une régression linéaire. En testant quels paramètres sont significativement non nuls, nous pouvons en déduire quelles variables explicatives influent sur l'expression du gène. Lorsque nous avons affaire aux expressions de nombreux gènes, nous sommes alors confrontés au problème des tests multiples : même en présence d'une seule variable explicative, il faut réaliser un test de nullité d'un paramètre pour chaque gène.

Au lieu de contrôler la probabilité d'avoir un faux positif, nous pouvons alors être tentés de contrôler l'espérance du rapport entre le nombre de faux positifs et le nombre de tests positifs, aussi appelé False Discovery Rate ou FDR. Il est facile d'obtenir une série de tests de nullité de certains paramètres basée sur le FDR, si chaque test réalisé pour un gène donné est indépendant des tests correspondant à chacun des autres gènes.

Dans notre cas, cela n'est pas vrai. Nous sommes donc amenés à introduire un modèle plus complexe, afin de retirer une structure des données transcriptomiques et ainsi produire des tests indépendants. Cette approche, développée par [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF], est appelée FAMT pour Factor Analysis for Multiple Testing. Elle est présentée dans la Section 4.2 dans laquelle l'implémentation d'outils informatiques pour faciliter l'utilisation de cette méthode est également détaillée.

FAMT se restreint à des relations linéaires entre expressions de gènes et variables explicatives. Pour modéliser des relations non-linéaires, il est possible d'utiliser la méthode SIR (pour Sliced Inverse Regression), introduite par [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], qui permet d'estimer B 1 -Introduction générale dans le modèle univarié suivant : y = f (x ′ B, ε), où f est une fonction quelconque de R K+1 dans R, x est une variable aléatoire de dimension p, y et ε sont des variables aléatoires univariées, x et ε sont indépendants et B est une matrice p × K de rang K, avec K ≤ p. Ainsi présenté, ce modèle ne permet pas de gérer des expressions de plusieurs gènes. Toutefois, si nous avons à notre disposition les expressions de q gènes, que l'on note y (1) , . . . , y (q) , nous pouvons le remplacer par le modèle suivant : 1) ), . . .

     y (1) = f 1 (x ′ B, ε (
y (q) = f q (x ′ B, ε (q) ),
(1.1) où f 1 , . . . , f q sont des fonctions quelconques de R K+1 dans R et ε (1) , . . . , ε (q) sont des variables aléatoires univariées indépendantes de x. On appelle espace EDR, l'espace engendré par B. Pour tout j = 1, . . . , q, en utilisant la méthode SIR, nous sommes en mesure d'obtenir un estimateur B j de B grâce à la j ème équation du modèle (1.1). Soit B = B 1 , . . . , B q , et soit Σ, l'estimateur usuel de la matrice de variance-covariance de x. Dans A new sliced inverse regression method for multivariate response regression, nous montrons que la matrice formée par les K vecteurs propres de la matrice B B ′ Σ correspondant aux K plus grandes valeurs propres converge en probabilité vers une matrice qui engendre l'espace EDR. Cet article constitue la Section 4.3 et illustre également le bon comportement de cette approche en simulation ainsi qu'avec deux jeux de données réelles.

Même lorsque l'on étudie les expressions d'un seul gène, la méthode SIR classique peut ne pas convenir à des jeux de données de trop grande taille. En effet, elle nécessite que le nombre p de variables explicatives à prendre en compte soit inférieur à la taille n de l'échantillon observé. Si cette condition n'est pas vérifiée, nous sommes alors dans un cas sous-déterminé. Le problème est que les estimateurs de B fournis par SIR sont généralement obtenus par une décomposition en éléments propres de la matrice Σ -1 M , où M est une matrice p × p calculée à partir des données et définie dans la Section 4.4. Or lorsque n < p, Σ n'est pas inversible.

Plusieurs approches ont été développées pour modifier la méthode SIR afin qu'elle puisse fonctionner dans ce cas. Par exemple, la méthode RSIR (pour Regularized SIR) de [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] est basée sur le remplacement de Σ par Σ(s) = Σ + sI p , où I p est la matrice identité de taille p×p et s est un réel positif. Le paramètre s optimal est déterminé par une technique de bootstrap. Nous proposons une méthode alternative appelée SIR-QZ qui choisit s comme étant le plus petit réel positif possible tel que l'algorithme QZ, qui calcule les éléments propres généralisés de M et Σ(s), renvoie un résultat dépourvu d'erreur numérique. Notez que la méthode SIR classique dépend d'un paramètre H : le nombre de tranches. Une autre fonctionnalité de SIR-QZ est de se servir de plusieurs estimations de x ′ B pour des valeurs de H différentes, afin d'obtenir une estimation qui dépend moins du choix de ce paramètre. Dans Comparison of sliced inverse regression approaches for underdetermined cases, RSIR et SIR-QZ sont comparés en simulation entre elles ainsi qu'avec deux autres approches. SIR-QZ y est également appliqué sur des données transcriptomiques. Cet article constitue la Section 4.4.

Chapter 2

General introduction

This thesis is made of modelling examples in biology, based on existing or original statistical tools and on datasets coming mainly from the EPOC laboratory (Environnements et Paléoenvironnements Océaniques et Continentaux) of the University of Bordeaux. After two introductory chapters, respectively in French and in English, Chapter 3 presents works about measurements of the activity of oysters. In Chapter 4, we describe statistical methods well-suited for transcriptomic data. We conclude this work in Chapter 5 where prospective research projects are also developed.

A stochastic model for bioindicators

In Chapter 3, we try to sum up the information which is contained in large samples. The signals we considered are indeed composed of high-frequency measures of the distance between both parts of the shell of an oyster. We call them valvometric signals. These data are for instance studied in order to give a diagnosis about the health of the animal. When several oysters from a single location are in poor health, changes in their environment are a plausible explanation. Finding a model to explain valvometric data is thus an important step toward the creation of efficient biosensors based on oysters to monitor the quality of water. We present at length the biological context related to these valvometric measurements in the introduction of As shown in Figure 2.1, two states can be observed in a valvometric signal. We name these states "open" and "closed". When an oyster is in good health, biologists observe that the transitions from a state to another are linked with the tide. This is not generally the case when the animal does not adapt well to its environment. Using a renewal process which can visit the states "open" and "closed" and an estimate of its survival function, we show in A hidden renewal model for monitoring aquatic systems biosensors (see the list of works) how we can perform a classification on some animals which only relies on their valvometric signal. This classification appears to be consistent with some information we have about their health. This article forms Section 3.1.

The renewal process is piecewise constant and the observed signal (Figure 2.1) is continuous. This signal also exhibit areas where it varies fast toward low values and quickly returns to a high level. These spikes are called microclosings and their minimum can sometimes be so low that it is difficult to separate microclosings from transitions from the open state to the closed state. In Section 3.1, we describe a method to estimate the renewal process which is hidden in the valvometric signal. This procedure allows us to determine a threshold. When the valvometric signal is below this threshold, we decide that the underlying process is in the closed state, otherwise it is in the open state. To find the threshold we need, we estimate a probability density function related to the valvometric signal. We assume that this density has two modes: one of them corresponds to the open state and the other one is linked to the closed state. We choose a density estimator which produces an estimate with two modes too. The location of the unique antimode of the estimate is then chosen to be the threshold.

Several methods were created to estimate a probability density function f when we have independent observations from the probability distribution defined by f at our disposal, and when we know how many modes f has. For example, we can apply a kernel density estimator or an estimator based on the level sets of f on these data. In Comparison of kernel density estimators with assumption on number of modes, we studied these estimators and another kernel density estimator which does not need any hypothesis about the number of modes of f , in a simulation study. This article forms Section 3.2.

The number of microclosings and their magnitude can also provide information about the health of an animal. It is however difficult to detect them because the valvometric signal can contain noise. In Section 3.3, we present an existing denoising method based on Haar wavelets. It was chosen among many others thanks to the simulation study from [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF]. We also studied in this section how to locate microclosings directly from the signal with noise.

To make these approaches easier to use, an implementation was realized with R and a user interface was created using the GTK+ library. These programs include additional features such as the ability to download valvometric data from a remote SQL server where they are stored as it is explained in Section 3.4. Some ideas are also presented in this section to optimize computational time when working with R.

Statistical methods to analyse genetic data

The biological problems we tackle in Chapter 4 concern the links between transcriptomic data and some other variables. The latter ones are named explanatory variables. Transcriptomic data are measures of how much genetic markers are expressed. We would like to explain how the explanatory variables make these expressions change. This is a regression context and the large size of the datasets is characterized either by a great number of variables to explain, which means a lot of genes, or by many explanatory variables. We provide two examples of such problems in Section 4.1.

We can link the changes in the expression of a single gene with the explanatory variables using a linear regression. By testing which parameters are significantly non null, we can deduce which explanatory variables have an effect on the gene expression. When we deal with expressions of several genes, we are facing multiple tests: Even if we only have one explanatory variable, we have to test one parameter for nullity, for each gene.

Instead of controlling the probability to obtain a false positive, we could prefer to control the expected value of the ratio between the number of false positives and the number of positive tests, also known as False Discovery Rate or FDR. If each test which is realized for a given gene is independent from the tests related to every other genes, we are able to produce a series of tests for nullity based on the FDR.

For our transcriptomic data, this assumption is not true. We then have to introduce a more complex model, in order to withdraw the structure of these data and then produce independent tests. This approach was developed by [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF] and is called Factor Analysis for Multiple Testing or FAMT. We describe it in Section 4.2 where we also give details of the user-friendly tool we implement for this method.

FAMT can only detect linear links between gene expressions and explanatory variables. The Sliced Inverse Regression (SIR) method can however find non-linear relations. It was introduced by [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], and aims at estimating B in the following univariate model:

y = f (x ′ B, ε),
where f can be any function from R K+1 to R, x is a p-dimensional explanatory variable, y and ε are univariate random variables, x and ε are independent and B is a p × K matrix with full rank and which satisfies K ≤ p. As it is, this model can not handle expressions of several genes. Let y (1) , . . . , y (q) be the expressions of q genes. We can replace the previous univariate model by the following one      y (1) = f 1 (x ′ B, ε (1) ), . . . y (q) = f q (x ′ B, ε (q) ),

(2.1)

where f 1 , . . . , f q can be any function from R K+1 to R and ε (1) , . . . , ε (q) are univariate random variables which are independent from x. The vectorial space spanned by B is called EDR space. For all j = 1, . . . , q, we are then able to obtain an estimate B j of B, using the SIR method and the j th equation of model (2.1). Let B = B 1 , . . . , B q , and let Σ be the usual estimator of the variance-covariance matrix of x. In A new sliced inverse regression method for multivariate response regression, we show that the matrix made of the K eigenvectors of the matrix B B ′ Σ which correspond to the greatest K eigenvalues converges in probability toward a matrix which spans the EDR space. This article forms Section 4.3 and also illustrate the good behaviour of this approach with analyses on simulated data and on two real datasets.

Even when we study expressions from a single gene, the usual SIR method may not suit too large datasets. More precisely, it needs that the number p of explanatory variables to take into account is less than the sample size n. If this requirement is not met, we have to deal with an underdetermined case. The issue is that estimates of B are generally provided by an eigendecomposition of the matrix Σ -1 M , where M is a p × p matrix which is computed from the data and defined in Section 4.4. When n < p, Σ is however not invertible.

Several improvements of the SIR method where then developed to make it work in this case. For example, the RSIR method (for Regularized SIR) from [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] replaces Σ by Σ(s) = Σ + sI p , where I p is the p × p identity matrix and s is a positive real. The optimal value of s is determined with a bootstrap technique. We also propose an alternative approach, called SIR-QZ. It is based on the QZ algorithm which computes generalized eigenelements of M and Σ(s). It chooses s as small as possible such that the QZ algorithm returns a result without any numerical error. Note that the usual SIR method relies on a parameter H: the number of slices. Another feature of SIR-QZ is to take advantage of several estimates of x ′ B for various values of H in order to produce a final estimate which is less sensitive to the choice of H. In Comparison of sliced inverse regression approaches for underdetermined cases, RSIR and SIR-QZ are compared with each other and with two other methods in a simulation study. SIR-QZ is also applied to some transcriptomic data. This article forms Section 4.4.

Chapter 3 A stochastic model for bioindicators

We propose in this chapter statistical methods dedicated to the transformation of a crude signal from an animal into a conclusion about its environment. Here, the environment is an aquatic one, the animals are bivalves and the measurement technique is the valvometry. A bivalve is a mollusc with a shell made of two articulated parts named valves. Given a location on each of these valves, a valvometric signal is then a sequence of measurements of the distance between them. The MolluSCAN Eye team currently composed of the engineers and researchers Pierre Ciret, Gilles Durrieu, Jean-Charles Massabuau, Mohamedou Sow and Damien Tran, studied such signals since 1997. This allows them to obtain insights about bivalve's physiology and aquatic toxicology. This team is supported by the EPOC laboratory (UMR CNRS 5805).

While in the following sections we will only focus on the measures they took from oysters, they have installed several valvometers around the world with various species. The signals we will consider come from Locmariaquer, Arcachon and Santander. The first two places are in France. The last one is in Spain. Another device was disposed in New Caledonia and worked from 2007 to 2009. It contained giant clams. Scallops were placed in the Svalbard archipelago and in Russia, where a mussel-based valvometer is also present. Valvometers rely on electrodes positioned on each valve of the animal as in Figure 3.1. More technical details about these devices can be found in Section 3.1.2 or at http://molluscan-eye.epoc.u-bordeaux1.fr/.

In Section 3.1, we present a suitable model for the valvometric data in order to diagnose the health status of the animals. This model relies on a kernel density estimator and on an assumption about the number of modes of a probability density function. We compare this kind of estimator with other ones in Section 3.2 in order to find a precise estimation of the location of the lowest local minimum of the probability density function.

Sections 3.1-3.2 aim at finding an overall diagnosis about the health status of oysters. Other interesting works concern more precise indicators. For example, [START_REF] Schwartzmann | In situ giant clam growth rate behavior in relation to temperature: A one-year coupled study of high-frequency noninvasive valvometry and sclerochronology[END_REF] observe the growth of bivalves using valvometers and another method called sclerochronology. [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica I. spawning reactions of the female and male[END_REF] studied spawning movements with an early type of valvometer. Scientists from EPOC witnessed the same behaviour with their installations. Spikes with a great amplitude are interesting features because bivalves only exhibit a lot of them in particular situations such as spawning or illness. They are for instance considered in [START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF] and in Section 3.3. In this section, we also tackle the issue of a possible noise in the signal. The analyses that lead to the model of Sections 3.1-3.2 needed a programming effort which is explained in Section 3.4.

A hidden renewal model for monitoring aquatic systems biosensors

This section is an article under revision and was written with Romain Azaïs and Gilles Durrieu. It was submitted in Environmetrics.

Introduction

Protection of the aquatic environment is a top priority for marine managers, policy makers, and the general public. Due to an increasing interest in the health of aquatic systems, there is a compelling need for the use of remote online sensors to instantly and widely distribute information on a daily basis. Among these sensors, bio-indicators are increasingly used and are highly effective to reveal the presence of low concentrations of contaminants through accumulation in tissues (see [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF][START_REF] Tran | Copper detection in the Asiatic clam Corbicula fluminea: Optimum valve closure response[END_REF]Tran et al., , 2007)). The ability of mollusk bivalves to adapt with the environment is one of the possible ways to assess water quality.

Monitoring their opening and closing activities over time is yet another method to evaluate the behavior of the bivalves in reaction to their environmental exposure.

The interest in investigating the bivalve's activities by recording their valve movements (valvometry) has been explored in ecotoxicology for more than 20 years. The basic idea of valvometry is to use the bivalve's ability to close its shell when exposed to a contaminant as an alarm signal [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF]Nagai et al., 2006 among others). Nowadays, valvometers are available on the market and use the principle of electromagnetic induction [START_REF] Sloff | Detection limits of a biological monitoring system for chemical water pollution based on mussel activity[END_REF][START_REF] Jenner | A new system for the detection of valve movement response of bivalves[END_REF] such as the Mossel Monitor [START_REF] Kramer | The valve movement response of mussels: a tool in biological monitoring[END_REF] or the Dreissena Monitor [START_REF] Borcherding | The dreissena monitor -first results on the application of this biological early warning system in the continuous monitoring of water quality[END_REF]. There has been a clear research interest in the recent years to measuring the bivalve's behaviors directly in real conditions [START_REF] Robson | Mussels flexing their muscles: a new method for quantifying bivalve behavior[END_REF][START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF][START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF].

These noninvasive valvometric techniques produce high-frequency data and different statistical models were built to analyze them [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF][START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF][START_REF] Jou | A dynamic artificial clam (corbicula fluminea) allows parcimony on-line measurement of waterborne metals[END_REF][START_REF] Jou | A dynamic artificial clam (corbicula fluminea) allows parcimony on-line measurement of waterborne metals[END_REF]. Here we propose a statistical procedure based on a four-state stochastic process to give inferences about oysters' health and provide some arguments about the healthiness of their environment. This method also exhibits a link between the tide and oysters' behavior, as shown in [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF].

The paper is organized as follows. In Section 3.1.2, we give a detailed presentation of the real data problem. Briefly, we consider valvometric data samples collected by the laboratory Environnements et Paléoenvironnements Océaniques et Continentaux (EPOC, http://www.epoc.u-bordeaux.fr/). In Section 3.1.3, we describe the state either open, or close, or in transition from open to close, or close to open, of an oyster by a fourstate renewal model. For this class of continuous-time stochastic processes, the path is piecewise-constant over time and changes at random times. In our particular framework, the motion is characterized by the conditional jump rate λ or by its integrate version Λ, called the cumulative jump rate. We provide a method to estimate Λ from a single trajectory of the renewal process. However, this process is not directly observed, but instead we have a signal that relies on it. Section 3.1.4 is devoted to extract from this signal the period of time when the oyster is in each state of the process. For this purpose, we estimate a density related to the distance between both parts of its shell, each day. This density is reasonably assumed to have two modes, the first one corresponds to a closed state of the shell and the second one to an open state. It leads to a threshold that enables us to approximate the time when oysters are in closed or open state. Thus, Section 3.1.3, together with Section 3.1.4, provides estimates of the survival functions related to both states. In Section 3.1.5, we apply the results of previous sections to biosensors for environmental monitoring. We especially exhibit links between groups of oysters and features related to the survey of the environment in different experimental sites, such as temperature variations and knowledge about their imminent death. Concluding remarks are given in Section 3.1.6.

Biological data

We first describe the experimental site and the animal species. Then we give some details on evaluation of valve activity. Afterwards, we provide information on data collection and transmission. Finally, exploring typical features of these environmental data, we explain which inferences are valuable from the biological point of view.

Data acquisition

The three monitoring sites we considered are respectively located in France, either in the Bay of Arcachon at the Eyrac pier (Latitude: 44 • 40 N, Longitude: 1 • 10 W), or at Locmariaquer (Latitude: 47 • 57 N, Longitude: 2 • 94 W) or in Spain, at Santander (Latitude: 43 • 44 N, Longitude: 3 • 79 W). A group of sixteen Pacific oysters, Crassostrea gigas, measuring from 8 to 10 cm length, are installed on each site. Every oyster has almost the same age (1.5 years old). They were placed in a traditional oyster farmer bag.

The electronic principle of monitoring was described by [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF] and further modified by [START_REF] Chambon | Influence of the parasite worm Polydora sp. on the behaviour of the oyster crassostrea gigas: a study of the respiratory impact and associated oxidative stress[END_REF] to handle difficult environmental conditions such as strong water current conditions and even storms. Some information about these specific aspects can be found on http://molluscan-eye.epoc.u-bordeaux1.fr. The main challenge was to ensure the complete autonomy of the system without in-situ human intervention for at least one full year. In brief, each animal is equipped with two light coils (sensors), of approximately 53 mg each (unembedded), fixed on the edge of each valve. These coils measure 2.5×2.5×2 mm and were coated with a resin sealing before fixation on the valves. One of the coils emits a high-frequency, sinusoidal signal which is received by the other coil. In each group of sixteen animals, one measurement is received every 0.1s (10 Hz). This means that each animal's behavior is measured every 1.6s. Every day, a data set with 864,000 pairs of values (1 distance value, 1 stamped time value) is generated. A first electronic card manages the electrodes and is in a waterproof case next to the animals. A second electronic card handling the data acquisition and the programmed emission is also in the field but outside the water (on a pier, a buoy, or in a light house). This unit is equipped with a GSM/GPRS modem and uses Linux operating system for driving the first control unit, managing the data storage, accessing the Internet, and transferring the data.

A self-developed software module runs on mobile phone technology. After each 24h period (or any other programmed period of time), the data collected are transmitted to a remote central workstation server where they are stored in text files. Every day, files from each site are inserted in a SQL database. This database is accessible with the software R

The biological issue

These measurements produce some characteristic features that can be examined in Figure 2.1. Oysters generally follow a two-stated pattern of behavior. The first state occurs when the distance between the parts of their shell is minimal and roughly constant, and the other one when this distance is bigger, varies more and is scattered with fast and short changes called microclosings. To date, link between tide and state of the oysters was shown [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF][START_REF] Tran | Field chronobiology of a molluscan bivalve: How the moon and the sun cycles interact to drive oyster activity rhythms[END_REF] in the sense that they are in open state when the tide is high and in closed state when it is low. However, oysters can close its shell for other reasons such as defense against an inhospitable environment. For example, pollution can affect oysters (see for instance [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF][START_REF] Adler | R in a Nutshell[END_REF]and Sow et al., 2011). Thus, detecting these changes of state can provide insights about the health of oysters. Because bivalves that we focus on are penned in where tide is semidiurnal, it makes sense to obtain at least one closed state and one open state every day like in Figure 2.1.

Four-state renewal model

In this section, we propose to model the behavior of oysters over time by a marked continuous-time process with four states. Indeed as shown in Figure 3.2, at each time one may consider that an oyster is either open, or close, or in a transition state from open to close or close to open. Among these four states, the open one and the closed one are of particular interest, given the remarks of Section 3.1.2.2. We thus consider the state space given by E = {O, T OC , T CO , C}, where O stands for open, while C stands for close and T OC (respectively T CO ) models the transition state from open to close (respectively from close to open). Finally, at each time t, X t ∈ E models the state of an oyster.

The process (X t ) is assumed to be a marked renewal process on E. This is a class of piecewise-constant stochastic models which change their location at random times, called jump times. Here, (T k ) denotes the sequence of the jump times of (X t ), where T 0 < 0 is unknown. For renewal processes, one often considers the inter-jumping times S k 's given, for any k ≥ 1, by S k = T k -T k-1 . As a consequence, the first inter-jumping time S 1 is unknown. One may also take into consideration the sequence (Z k ) of the locations of (X t ). For any integer k, Z k = X T k .

The relation between the continuous-time process (X t ) and the discrete-time one (Z k , T k ) is given by

X t = Z k for T k ≤ t < T k+1 .
As a consequence, the discrete-time process (Z k , T k ), and equivalently the Markov chain (Z k , S k ), contains all the information of (X t ). In our particular four-state case, the motion of the process (X t ) depends only on its conditional jump rate, which determines the distribution of the S k 's. The dynamic of (Z k , S k ) is defined, for all k ≥ 0, by

• if Z k = O, Z k+1 = T OC , • if Z k = T OC , Z k+1 = C, • if Z k = C, Z k+1 = T CO ,
• and if

Z k = T CO , Z k+1 = O,
and for all k ≥ 1 and t ≥ 0,

P(S k+1 > t|Z k , . . . , Z 0 , S k , . . . , S 1 ) = P(S k+1 > t|Z k ) = exp - t 0 λ(Z k , s)ds .
An example of a path of (X t ) is given in Figure 3.2. The function λ : E × R + → R + is called the (conditional) jump rate of the process (X t ). One may consider also the cumulative jump rate Λ :

E × R + → R + which is the integrate version of λ, ∀(x, t) ∈ E × R + , Λ(x, t) = t 0 λ(x, s)ds.
Each of both these functions fully characterizes the conditional distribution of the interjumping times.

Here, we investigate the nonparametric estimation of the cumulative jump rate for an oyster from the observation of its behavior. We consider a nonparametric framework because no prior information about the time spent by an oyster in a given state allows us to choose a parametric family of distributions. We shall obtain for each considered oyster a function which completely characterizes the distribution of the period spent in the state open or in the state close. In the sequel, m denotes the number of observed jumps. The observation of the continuous-time process within one week corresponds with the observation of about eighty jumps.

For any t ≥ 0, we propose to estimate the cumulative jump rate Λ(x, t), x ∈ E, by the Nelson-Aalen estimator Λ m (x, t) (see Andersen et al., 1993 and the references therein) given by

Λ m (x, t) = m-1 k=1 R m (x, S k+1 )1 {Z k =x} 1 {S k+1 ≤t} , (3.1)
where 1 A takes the values 1 and 0 according to whether the condition A is satisfied or not and

R m (x, t) = 1 Lm(x,t) if L m (x, t) > 0, 0 otherwise, with L m (x, t) = m-1 k=1 1 {Z k =x} 1 {S k+1 ≥t} .
We exclude the first inter-jumping time S 1 from the estimation procedure since it is unknown.

The Nelson-Aalen estimator is well defined in the framework of the multiplicative intensity model, introduced by Aalen in the seventies. Fortunately, Azaïs et al. have shown that this model is well suited for the particular case of marked renewal processes with discrete state space. In this context, we establish in Proposition 1 the consistency of the estimator under some additional conditions imposed on the asymptotic behavior of R m and L m when m goes to infinity. We also prove in Lemma 2 that these conditions are verified under a reasonable assumption on λ.

Proposition 1. Let (x, t) ∈ E × R + . Assume that the conditions ∀ε > 0, lim m→∞ P t 0 R m (x, s)λ(x, s)ds > ε = 0, (3.2)
and ∀ε > 0, lim

m→∞ P t 0 1 {Lm(x,s)=0} λ(x, s)ds > ε = 0, (3.3)
are verified. Then, we have When the process (X t ) is hidden, one may only obtain some approximations S k 's of the inter-jumping times S k 's. For every k ∈ {2, . . . , m}, let

∀ε > 0, lim m→∞ P sup 0≤s≤t Λ m (x, s) -Λ(x, s) > ε = 0. ( 3
S k ∈ [S k , S k + ζ],
for a small constant ζ > 0 . We do not compute the Nelson-Aalen estimator Λ m (x, t) but an approximation of this estimator Λ m (x, t) from the S k 's. Since, for all the integers k,

R m (x, S k+1 ) = R m x, S k+1 , we have ∀x ∈ {O, C}, ∀t ∈ R + \ k S k -ζ, S k , Λ m (x, t) = Λ m (x, t).
Thus, the error between the Nelson-Aalen estimator and its approximation is equal to 0 for most t ≥ 0.

Besides the cumulative jump rate, one can also estimate the conditional survival functions H associated with Λ. Indeed, these functions are easier to classify than the Nelson-Aalen estimators because they take values between 0 and 1, whereas the range of values taken by Λ m depends on m. The Fleming-Harrington estimator [START_REF] Fleming | Nonparametric estimation of the survival distribution in censored data[END_REF] of H can easily be computed from our modified Nelson-Aalen estimator. For any x ∈ E, for any t ≥ 0, it is given by

H m (x, t) = exp -Λ m (x, t) .
By considering the process t → 1 {Xt=T OC } , one may also estimate the survival function of the time of the first return in the state T OC . Let J m be the subset of {2, . . . , m -4} such that for any k ∈ J m , Z k = C. Using, for each k ∈ J m , an approximation S k+1 of S k+1 +S k+2 +S k+3 , one may simply compute an approximation of the Fleming-Harrington estimator of the survival function of the time of the first return in the state T OC . This estimator is denoted H T OC m .

Estimation procedure

In the application, since we do not observe the process (X t ), the estimator Λ m (x, t) of Λ(x, t) can not be computed. We assume that we observe instead a stochastic process (G t ) that looks like the signal in Figure 2.1. This process should behave like (X t ) by taking low values when X t = C, high values when X t = O and intermediate values otherwise. The aim of this section is to use (G t ) to construct estimators S k of the observations of S k for some appropriate values of k. We would then be able to compute Λ m , H m and H T OC m . In Section 3.1.4.1, we link (G t ) with a probability density function f that is designed to summarize the information contained in (G t ). We especially assume that f has a given number of modes. In Section 3.1.4.2, we propose a modification of kernel-based estimators of f for equally spaced sample.

A probability density for oysters' openings

For a fixed time interval and a fixed oyster, consider that (G t ) is the process that returns the oyster's opening amplitude for a given instant. The process (G t ) is studied for t ∈ [0, T max ] where T max is a positive real. It takes values in an interval [a, d] ⊂ R + . For all ω in the probability space Ω, let G be a function from [0, T max ] to [a, d] defined by

G : t → G(t) = G t (ω)
. The assumptions about G made in this article are thus made for all ω ∈ Ω. We assume that G is continuously differentiable except for a countable set of points on which G is only continous. For all i ∈ {1, . . . , n}, let V i be a random variable that follows the distribution U [0,Tmax] and Y i the random variable defined by Y i = G(V i ). Let f denote the continuous probability density of each Y i given G. Such a density can be also defined as the occupation density of G on [0, T max ]. See Geman and Horowitz (1980, § 2) and Bosq and Blanke (2008, Section 9.2) for details.

The following assumptions allow us to link f and G(t) for t such that X t ∈ {T OC , T CO }.

Assumptions 3.

(i) It exists a couple (b, c), with a < b < c < d, such that,

∀t ∈ [0, T max [, X t = C ⇒ G(t) < b and X t = O ⇒ G(t) > c. (ii) inf {|G ′ (t)| : G(t) ∈ [b, c]} > 0 and G ′ is continuous on {t : G(t) ∈ [b, c]}.
Assumptions 3 are reasonable with respect to the signal observed in Figure 2

∀p ∈ [b, c], f (p) = 1 T max k∈Im(p) 1 |G ′ (u k )| , (3.5)
where m = max{k : T k < T max } and I m (p) is the subset of {1, . . . , m} made of integers k such that there exists a unique

u k ∈ [T k , T k+1 [ which satisfies G(u k ) = p. For every x ∈ R, if in points t ∈ [0, T max ]
where G(t) = x, the values of |G ′ | are overall great, f (t) is low and conversely. Let N (f ) be the number of modes of the density f . We postulate that the values of |G ′ | are lower when the oyster is open or closed than when it is in a transition state. We thus assume that N (f ) = 2 which is reasonable for a period of time [0, T max ] throughout which the oyster visits both open and closed state. We also assume that during transition states |G ′ (t)| is greater for t such that G(t) ∈ [b, c] than for t such that G(t) is close to b or c but not in [b, c]. This remarks leads to Assumptions 4.

Assumptions 4.

(i) m ≥ 3.

(ii) N (f ) = 2, f has no flat part and has a unique antimode located in θ.

(iii) ∀t ∈ {b, c}, ∀k ∈ I m (t), ∃ε k > 0, ∀ε ∈]0, ε k [, ∃v k ∈ {u k -ε, u k + ε}, G(v k ) / ∈ [b, c] and |G ′ (v k )| ∈ ]0, |G ′ (u k )|[ , where u k is the unique point in [T k , T k+1 [ such that G(u k ) = t.
Note that it is easy to take T max great enough so that the first point of Assumptions 4 is satisfied. Given Assumptions 3 and 4 and equation (3.5), we then have

P (θ ∈ [b, c]) = 1.
In this section, we showed that we can link the stochastic function G that describes an oyster's opening during a given time interval and a probability density f . In Section 3.1.4.2, we propose some estimators for f and θ. The estimate of θ will then be taken as a threshold and the instants when G crosses it will lead to approximations of the S k 's.

Density estimation procedure

For t ∈ R, we use the following kernel density estimator f n (t) of f (t)

f n (t) = 1 nh n n i=1 K t -G i Tmax n h n , (3.6)
where K is the Gaussian kernel defined, for every t ∈ R, by K(t) = 1 √ 2π e -1 2 t 2 and h n is a positive real parameter called the bandwidth. Some pointwise properties about kernel estimators for occupation densities are available in [START_REF] Wang | Asymptotic theory for local time density estimation and nonparametric cointegrating regression[END_REF].

The estimator f n allows us to estimate θ under the following assumptions on the bandwidth h n and G.

Assumptions 5.

(i) f is uniformly continuous.

(ii) h n → 0 and nh 2 n → ∞ as n goes to infinity.

(iii) G ′ is bounded.

(iv) f n has a unique antimode located in θ n .

Under the second point of Assumptions 4 and Assumptions 5, we have

P lim n→∞ θ n = θ = 1.
Note that because G is measured from a living being, the third point of Assumptions 5 is reasonably satisfied.

To ensure the convergence of θ n toward θ, the bandwidth h n remains to be chosen. To satisfy the last point of Assumptions 5, mimicking [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF], we choose h n = h crit , which is defined as

h crit = min{h : N ( f n ) ≤ 2}.
When K is the Gaussian kernel, [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF] showed that h crit can be easily computed.

For any t ∈ R, define f n as

f n (t) = 1 nh n n i=1 K t -Y i h n . (3.7)
The second point of Assumptions 5 for h crit is convincing because if K is the Gaussian kernel, we have the following points.

• Let h crit be defined as h crit using f n instead of f n . We have for any ε > 0 that P( h crit > ε) → 0 and that P(n h 2 crit > ε) → 1 when n goes to infinity (see Mammen et al., 1991 and Section 3.2).

• for such a bandwidth, f n converges uniformly in probability toward f (see Devroye and Wagner, 1980 and Section 3.2).

• For a given random bandwidth h rand used in both f n and f n and such that for any ε > 0, P(h rand > ε) → 0 and P(nh 2 rand > ε) → 1 when n goes to infinity, we can prove that f n converges uniformly in probability toward E f n .

When n is large enough, we are now able to build a threshold θ n from G which is only crossed by G at points t ∈ [T k , T k+1 [ such that Z k ∈ {T OC , T CO }. In addition, asymptotically, G crosses only once θ n in each of such intervals [T k , T k+1 [, because of the second point of Assumptions 3. For every

k ∈ I m ( θ n ) let t k be the unique point in [T k , T k+1 [ such that G(t k ) = θ n . Thus, k ∈ I m ( θ n )\{max{I m ( θ n )}}, implies that k + 2 ∈ I m ( θ n ) and we have t k+2 -t k ∈]S k+2 , S k+2 + ζ[, where ζ = 2 max k∈Im( θn) {S k+1 }. The value of ζ is small because when k ∈ I m ( θ n ), Z k ∈ {T OC , T CO }. We then set S k+2 = t k+2 -t k for all k ∈ I m ( θ n )\{max{I m ( θ n )}}.
Thus, for x ∈ {O, C}, and t ≥ 0, we are able to build the estimators Λ m (x, t), H m (x, t) introduced in Section 3.1.3. Note that ∀k ∈ J m , {k -1, k + 3} ⊂ I m ( θ n ) and we can define

S k+1 = t k+3 -t k-1 ∈]S k+1 + S k+2 + S k+3 , S k+1 + S k+2 + S k+3 + ζ[,
and then compute H T OC m (x, t).

Real data application

In this section, we apply results from the previous sections to oysters from different experimental sites. We give H m (x, t) for x ∈ {O, C} in Figure 3.3 for four oysters studied during 8 days. We computed a value of θ n for each day and each oyster and pooled the resulting values of S i in a single estimator Λ m (x, t). The first two animals, corresponding to solid and dashed lines do not have reason to be in poor health. The oyster from Santander (mixed line), was in a harbor in Spain which is an environment that might not be clean enough for this kind of animal, and the last one (dotted line) died the 20th of May 2009. In Figure 3.3(b), one can observe that estimated survival functions related to the oyster from Santander and to the dying one are close from each other. In addition they can be separated from the two others. The latter functions decrease less than the previous ones for small values of t. This means that the oyster of Santander and the dying one exhibit more short opening periods than the two others. A big amount of such periods is characteristic of oysters that live in an inhospitable environment. Groups of curves of Figure 3.3(b) can also be noticed in Figure 3.3(a), but this is less obvious. Since the oyster from Santander sometimes spent more than 40 hours closed, the mixed line is far from 0 for t = 12. Another feature that can arise for oysters with health problems is that they do not open or do not close themselves during a long time. In this case, we can have m = 1 and N (f ) = 1 when T max = 24, which provides incorrect estimates H m (x, t). Hopefully, because of spikes in open state (see Figure 2.1), and because of electronic noise in closed state, the inaccurate θ n leads in both cases to small values of S i . Thus, the corresponding function H m (x, t) is close to estimates for oysters that exhibit a lot of microclosings. Because we aim at classifying oysters according to their health status, this does not bother us and we use H m (x, t) even if m = 1 and N (f ) = 1. Remark 1. An alternative way to handle a hidden process with a finite number of states and an observed one taking real values is the Hidden Semi-Markov Models (HSMM). O' Connell and Højsgaard (2011) implement in the R package mhsmm a procedure to estimate X t for t ∈ [0, T max ] in this context. They assume that the distribution G t depends only on X t , which means that G t given X t is independent of G s for any s < t. In addition, while the distributions of the S k 's given Z k can be non-parametrically estimated, the distribution of G t given X t has to belong to a parametric family. For oysters' opening measured with high frequency, these assumptions are not reasonable and the method to approximate X t is time consuming. The corresponding algorithm takes more than 4 hours to handle the signal of an oyster for a given day while finding a threshold with our method requires a few seconds.

For each animal of Figure 3.3, the estimated survival function H T OC m of the time of the first return in T OC is plotted in Figure 3.4. Curves in solid line and in dashed line decrease fast when t is around 12h and we know that in the locations where the oysters are placed, the tide has a period of approximately 12.41 hours. Hence, oysters' openings that correspond to these curves seem to be tide-driven, as they should be when the animals are not perturbed (see [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF]. This let us think that these oysters are in good health.

It can be unclear why we take T max = 24 to compute θ n while we use data from 8 days for each oyster. Notice that the oysters grow and that the distance between the parts of their shell can be greater when they are closed than when they were open before. That is why N (f ) can not be assumed to be equal to 2 for a T max too large.

We studied 9 other oysters than the 4 previously considered, during 8 days. A value of θ n was computed for each oyster and for each day, which led to estimates H T OC m . Among these 9 animals, 2 of them died at most 2 weeks after the period when the data were collected, 5 of them had no known reason to be in poor health and 3 of them lived in the harbor of Santander when the measurements were taken. The last one, installed at Locmariaquer went through temperature variations that may have perturbed it. We set N = 100 and we measured H T OC m (t) for t ∈ 0, 1 N -1 T max , 2 N -1 T max , . . . , T max , for each oyster. We performed a principal components analysis (PCA) on the resulting N × 13 matrix. This first principal component explains a large part of the total variability of the data set (88%) and the second one explains only 9% of it. We then clustered the oysters into two groups by using a k-means algorithm on the scalar products between the values of H T OC m (t) of every oyster and the first principal component of the PCA. The scalar products corresponding to the first two components of the PCA are plotted in Figure 3.5. This graphic also presents the result of the k-means algorithm. In the second group, the dying oysters, those from Santander and the one that suffered temperature variations are gathered. These animals have reasons to be in poor health.

Concluding remarks

This paper presented a framework to analyze a renewal process with a finite number of states when this process is not directly observed. Instead, we have a signal that takes values in different intervals when the process is in different special states of interest. Theses intervals are assumed to be disjoint and the states that are not related to any interval are called transition states. The absolute value of the derivative of the signal should be large in these transition states but it is allowed to take great values in the other states too.

Focusing in a four-state model with two transition states and two states of interest, we described how such a renewal process can underly signals made of measures of the distance between the parts of the shell of oysters. We gave convergences of an estimator of the cumulative jump rate of an observed process. Then, we summarized the information contained in the signals with a bimodal probability density function that is estimated by a bimodal function. The location of the antimode of this estimate provide a threshold that separate out the instants when the oyster is in each states of interest. We thus are able to estimate the cumulative jump rate for the states of interest, from a signal of oysters' openings.

Finally, for a transition state, we studied an estimator of the survival function of the first return time that appears to be different for oysters that have a tide-driven behavior than for animals that have health problems or live in an inhospitable environment. A classification procedure on these functions succeeded in separating out both kinds of oysters. Because the quality of their environment and their health status are closely related, this approach seems to us to be a new step toward the use of valvometric devices as biosensors for water quality.

Comparison of kernel density estimators with assumption on number of modes

This section is an article which is accepted for publication in Communications in Statistics -Simulation and Computation. It was written with Gilles Durrieu and Jérôme Saracco.

Introduction

Since the seminal papers of [START_REF] Parzen | On estimation of probability density function and mode[END_REF] and [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF], the use of kernels to find an estimate fK,h of a density function f of a random variable X is widely studied because of the advantages of the nonparametric point of view. Let (X 1 , . . . X n ) be a vector of independent and identically distributed random variables generated from f . For t ∈ R, the kernel density estimator fK,h (t) of f (t) can be defined as

fK,h (t) = 1 nh n i=1 K t -X i h (3.8)
where K is the kernel and h is the bandwidth on which the amount of smoothness of fK,h relies. Most of the time, K is a probability density function and h is a positive real. The larger the bandwidth, the smoother the estimate is. The choice of the bandwidth h is an important area in kernel estimators research field. Even if there exists a sufficiently large interval around the optimal bandwidth where fK,h stays roughly the same [START_REF] Scott | Multivariate Density Estimation : Theory, Practice and Visualization[END_REF], p. 161), h needs to be carefully determined. To perform this choice, one can use biased or unbiased cross-validation [START_REF] Rudemo | Empirical choice of histograms and kernel density estimators[END_REF]; [START_REF] Scott | Biased and unbiased cross-validation in density estimation[END_REF]) as well as plug-in methods [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF]) among other approaches.

In practical situations, the scientist that brings the data to analyze is able to determine if the estimated density function is smooth enough or not. In this paper, we are interested in using this information on the necessary amount of smoothing in order to set the corresponding bandwidth h for the estimator fK,h . More precisely, we will assume a fixed number N (f ) of modes of f . We will introduce and study the bandwidth h crit,k which is the smallest one such that the estimator fK,h has at most k modes. Thus, the definition of h crit,k is:

h crit,k = min N ( fK,h )≤k h, (3.9)
for a given k ∈ N * . The minimum is well-defined for the different kernels we will consider.

The link between h and N ( fK,h ) has been studied by several authors. With a Gaussian kernel (i.e. K is the density function of the standardized normal distribution), according to [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF], the function h → N ( fK,h ) is decreasing, which allows him and [START_REF] Mammen | Some asymptotics for multimodality tests based on kernel density estimates[END_REF] to test the number of modes of f . For many other kernels among those with bounded support, we do not have these kind of properties, but we have at our disposal a visualization tool called the "mode tree" (see for details [START_REF] Minnotte | The mode tree: A tool for visualization of nonparametric density features[END_REF] or [START_REF] Minnotte | The Bumpy Road to the Mode Forest[END_REF]). Other theoretical results are also available in the literature, see for instance [START_REF] Hall | Bump hunting with non-gaussian kernels[END_REF]. A method to find an estimate of f based on this kind of assumptions already exists (Polonik (1995a)).

The paper is organized as follows. In Section 3.2.2, we give asymptotic results for the density kernel estimator fK,h crit,k where K is the Gaussian kernel. We also present theoretical results for the uniform kernel. Define h crit = h crit,N (f ) . In Section 3.2.3, we present a simulation study which compares fK,h crit to another kernel estimator that is taken as a benchmark and relies on a widely-used bandwidth. The procedure from Polonik (1995a) was also considered in simulation to try to determine the better approach to include information about the number of modes of f in its estimation. Then, in Section 3.2.4, we describe constraints under which one can use h crit,k in the context of mixture models. We apply it to environment monitoring data in Section 3.2.5. Lastly, concluding remarks are given in Section 3.2.6.

Estimating a density with N (f ) modes

In this section we study the kernel density estimator fK,h given in (3.8) with the bandwidth h crit,k defined in (3.9). For our purpose, we only consider two kernels:

• the uniform kernel defined for t ∈ R as K(t) = 1 [-1 2 , 1 2 ] (t),
• and the Gaussian kernel defined for t ∈ R as

K(t) = 1 √ 2π e -1 2 t 2 .
The first kernel has a bounded support, this is not the case for the Gaussian kernel. We also describe two alternatives of fK,h crit,k which are respectively fK,h SJ , where h SJ is the bandwidth given by [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] plug-in method (see Section 3.2.2.4) and Polonik (1995a) estimator based on density contour clusters (see Section 3.2.2.5).

3.2.2.1

Assumptions on the density f of X

We need the following assumptions on the density f of X in order to have a density with N (f ) modes which can be properly estimated.

(H1) f is uniformly continuous on R.

(H2) ∃(r, s), -∞ < r < s < +∞, f (x) = 0 ⇒ x ∈ [r, s] and x ∈]r, s[⇒ f (x) = 0. (H3) ∃! (z 1 , z 2 , . . . , z 2N (f )-1 ) ∈ ]r, s[ 2N (f )-1 , ∀i ∈ {1, 2, . . . , 2N (f ) -1}, f (1) (z i ) = 0 and sign(f (2) (z i )) = (-1) i ,
where f (q) is the q th derivative of f .

(H4) f ∈ C 2 (]r, s[).
(H5) lim t↓r f (1) (t) > 0 and lim t↑s f (1) (t) < 0.

(H6) ∀x ∈]r, s[ that satisfies f (1) (x) = 0, f (2) (x) = 0.
Remark 2. (H1) follows [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF]. The authors give an asymptotic result with the L ∞ norm, that we discuss in Section 3.2.2.3. (H2) -(H6) are taken from [START_REF] Mammen | Some asymptotics for multimodality tests based on kernel density estimates[END_REF].

A computable bandwidth

For the Gaussian kernel, some interesting results on h crit,k already exist. They underline that the bandwidth h crit,k is easily computable. Indeed, [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF] shows that the function h → N ( fK,h ) is decreasing and right continuous. This ensures computability of h crit,k with the desired accuracy by a dichotomous search. With the assumption that

h crit,k ∈ [h 1 , h 2 ]
, and if we want to obtain it with an error less than h 2 -h 1 2 ň , we have to compute N ( fK,h ), ň times, for various h. If for each h, to determine N ( fK,h ), fK,h is computed in ñ points, then the computational complexity of the whole algorithm to find h crit,k is equal to O(nñň), where n is the sample size. In our simulations, we often take ñ = 10000 and ň ≤ 30.

For the uniform kernel, we provide a similar result by explaining that h → N ( fK,h ) is piecewise constant and has at most n(n-1) 2 jumps. Besides, if X (i) i∈{1,...,n} is a set of ordered random variables from which we want to compute fK,h crit,k , the locations of the jumps are equal to X (j) -X (i) for some i ∈ {1, . . . , n -1} and some j ∈ {i + 1, . . . , n}. This means that we are able to find h crit,k by analyzing values of fK,h between jumps.

Remark 3. The number of jumps in h → N ( fK,h ) is not bounded by n(n-1) 2 for every kernel. Indeed [START_REF] Hall | Bump hunting with non-gaussian kernels[END_REF] studied the set of points X(ω) = (-1, 0, 1) and drew

N ( fK θ ,h ), with K θ (x) = C θ (1 -x 2 ) θ 1 [-1,1] (x)
as a function of h and θ, where C θ ensures that K θ L 1 = 1. For example for θ = 1.5, one can find 4 different values in h → N ( fK θ ,h ), which is greater than n(n-1) 2 = 3.

Asymptotic results on fK,h crit,k

Proof of consistency for this estimator toward f is not trivial since h crit,k is data-driven. However, for the Gaussian kernel, the pointwise convergence in probability was proved by [START_REF] Futschik | On the consistency of kernel density estimates under modality constraints[END_REF], when k ≥ N (f ). We have furthermore the uniform convergence of the estimate toward the true density in probability and the convergence of the integrated absolute error to 0 in probability. To explain these results, we first find conditions for a given data-driven bandwidth h n under which some asymptotic properties can be shown for fK,hn . This is realized in Theorem A by combining Theorem 2 of [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF] about the L ∞ distance between fK,hn and f and Theorem 3.3 from Devroye (1987, p. 38), concerning the L 1 distance.

Theorem A [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF], [START_REF] Devroye | A Course in Density Estimation[END_REF]). Let f be a probability density satisfying (H1), h n a random bandwidth depending on X. If we assume the following hypotheses:

(F1) K is a Riemann integrable probability density, (F2) sup x∈R K(x) < ∞, (F3) ∞ 0 sup |x|≥z K(x)dz < ∞, (F4) ∀ε > 0, P(h n > ε) → 0, when n → ∞, (F5) ∀A > 0, P(nh n > A) → 1, when n → ∞, then, we have, for n → ∞, P sup t∈R fK,hn (t) -f (t) > ε → 0, and P R fK,hn (t) -f (t) dt > ε → 0.
The following theorem states that Theorem A can be applied with h crit,k , when K is the Gaussian kernel and k ≥ N (f ).

Theorem 6. Let f be a density satisfying (H1) -(H6) and let fK,h crit,k be the estimator of f with the Gaussian kernel K and the bandwidth h crit,k given in (3.9), with k ≥ N (f ). Then we have, for n → ∞,

P sup t∈R fK,h crit,k (t) -f (t) > ε → 0, and P R fK,h crit,k (t) -f (t) dt > ε → 0.
Proof. When K is the Gaussian kernel, (F1) and (F2) are satisfied as well as (F3) because K is defined and decreasing on [0, ∞[ and then

∞ 0 sup |x|≥z K(x)dz = ∞ 0 K(x)dx < ∞.
For this kernel, Corollary 2.1 from [START_REF] Mammen | Some asymptotics for multimodality tests based on kernel density estimates[END_REF]

implies that if k ≥ N (f ), P n -1/4 ≤ h crit,k ≤ n -1/6 → 1,
and we thus have (F4) and (F5) (see also Lemma 1 and Lemma 2 from [START_REF] Futschik | On the consistency of kernel density estimates under modality constraints[END_REF]). Theorem A can then be applied.

Assuming some regularity conditions on the kernel, [START_REF] Hall | Bump hunting with non-gaussian kernels[END_REF] proved similar results in their Theorems 3.1 and3.2 (pp. 2130-2131). These conditions on the kernel are stronger than continuity on R and thus the uniform kernel does not satisfy them. In that case, we prove that we cannot have (F4) in the following theorem.

Theorem 7. For any probability density function f of X, let fK,h crit,k be the estimator of f when K is the uniform kernel with h crit,k given in (3.9). Then we have h crit,k increasing with n, for all k ∈ N.

The proof is given in A.3.9. See also Appendix A.1 for details.

Sheather and Jones' plug-in method to choose a bandwidth

In the more general context of estimating a density without assumption on the number of its modes, algorithms that provide a suitable bandwidth h are of particular interest. Among a large selection of procedures, we focus on the plug-in method developed by [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] (see also [START_REF] Jones | Using non-stochastic terms to advantage in kernel-based estimation of integrated squared density derivatives[END_REF]), which leads to the bandwidth h SJ . We chose it because h SJ has good asymptotic properties and is easy to compute. This bandwidth is designed to minimize the asymptotic mean integrated squared error (AM ISE) between fK,h and f , defined as:

AM ISE(h) = R(K) nh + 1 4 σ 4 K h 4 R(f (2)
), (3.10) where for any function ψ, R(ψ) = ∞ -∞ ψ(x) 2 dx and σ 2 K is the variance of a random variable of density K. To minimize AM ISE(h), R(f (2) ) must be estimated. For this purpose, modifying an estimator studied by [START_REF] Hall | Estimation of integrated squared density derivatives[END_REF], [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] used the following one:

R(f (2) ) = 1 n(n -1) h5 n i=1 n j=1 K(4) X i -X j h , (3.11)
where K is allowed to be different from the kernel K used in the estimate of f , K(4) is the fourth derivative of K and h = 2K (4) (0) 3) ) is similar to the one in (3.11). It requires a new bandwidth ȟ chosen to be equal to 0.912 λn -1/9 , where λ is the sample interquartile range. Finally, [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] showed that

nσ 2 K R(f (3) ) 1/7 . Estimator R(f (3) ) of R(f (
h SJ = arg min h R(K) nh + 1 4 σ 4 K h 4 R(f (2) ) . Let h opt = arg min h (M ISE(h)),
h SJ h opt = 1 + O P (n -5/14 ).
This means that h SJ is close to the bandwidth that minimizes the expected L 2 distance between a kernel density estimator and the true density. Note that, for the Gaussian kernel, this result is valid for a density with three derivatives and which satisfies for all x and y:

∃Z > 0, f (3) (x) -f (3) (y) ≤ Z |x -y| 1 4 .
Remark 4. Among other methods to estimate R(f ) in (3.10), one can replace f by the density of a normal distribution as described in Silverman (1986, Section 3.4.2). This is called the normal reference rule. Another alternative is to use the density of a mixture of Gaussian distributions instead of f , in R(f ). This technique was introduced by [START_REF] Einbeck | A number-of-modes reference rule for density estimation under multimodality[END_REF]. It relies on several approximations but the resulting bandwidth is fast to compute and provides small values of M ISE(h) in simulation. The authors explain that the number m of components of the underlying Gaussian mixture should be set with respect to N (f ). Notice that the number of modes of the estimate which is produced by this procedure is not necessarily equal to m.

Polonik's estimator based on excess mass location

The kernel density estimator fK,h aims to associate a fixed point t with a value fK,h (t) as close as possible to f (t). Another approach, described in this subsection, tries to determine for every given λ ∈ [0, ∞[, the set Γn,C (λ) which is the most similar to Γ(λ) = {t : f (t) ≥ λ}, where C is a set of unions of disjoint intervals of R chosen such that for every λ ∈ [0, ∞[, Γ(λ) lies in C. In our case, C is made of every unions of at most N (f ) disjoint intervals. This procedure was developed by Polonik (1995b,a) and is related to the test for multimodality introduced by [START_REF] Müller | Excess mass estimates and tests for multimodality[END_REF]. It leads to an estimator of

f (t) = ∞ 0 1 Γ(λ) (t)dλ defined by: fP (t) = ∞ 0 1 Γn,C (λ) (t)dλ,
where Γn,C (λ), the so-called empirical generalized λ-cluster in C, satisfies Γn,C (λ) = arg max

{C(i)} i∈{0,...n} i n -λµ(C(i)) .
For every Polonik (1995a, Theorem 3.1) gives that:

C ∈ C, µ(C) is the sum of the length of all disjoint closed intervals in C. C(i) is one of the narrowest elements of C(i), which means that ∀C ∈ C(i), µ(C(i)) ≤ µ(C), where C(i) is the subset of C such that ∀C ∈ C(i), n j=1 1 C (X j ) = i. Then, if R fP (t)dt = 1,
P lim n→∞ R fP (t) -f (t) dt = 0 = 1.
We have R fP (t)dt = n-1 n , but, the proof of Theorem 3.1 from Polonik (1995a) still works for R fP (t)dt ≤ 1. Thus, to obtain a value for fP (t), C(i) should be determined for every i ∈ {0, . . . , n}. Then, because Γn,C (λ) is the C(i) that maximizes n + 1 linear functions, it only changes for a finite number k n of λ, with k n ≤ n. Let λ (i) i∈{1,...,kn} be the set of this ordered change points and take λ (0) = 0. When λ ≥ λ (kn) , Γn,C (λ) = C(0) and µ(C(0)) = 0, and we can set Γn,C (λ) = C(0) = ∅. Hence, another expression for fP is given by

∀t ∈ R, fP (t) = kn-1 i=0 λ (i+1) -λ (i) 1 Γn,C (λ (i) ) (t) ,
which can be used to computed fP (t).

Simulation study

In this simulation study, we compare four density estimators: fK,h crit based on the bandwidth h crit with both Gaussian and uniform kernels, fK,h SJ based on the Sheather and Jones' bandwidth with the Gaussian kernel, and Polonik's estimator fP . We specifically show numerical convergences of fK,h crit with the Gaussian kernel and illustrate consequences of the unsatisfied requirement of the uniform kernel.

Simulated data and quality assessment of the estimates

To generate simulated datasets, we use a beta mixture, a Gaussian mixture and a mixture of t-distributions. The beta mixture model is defined by:

X ∼ B(α 1 , β 1 ) with probability p 1 , B(α 2 , β 2 ) with probability p 2 = 1 -p 1 . (3.12)
Note that the corresponding density of X is:

f 1 (t) = p 1 Γ(α 1 + β 1 ) Γ(α 1 )Γ(β 1 ) t α 1 -1 (1 -t) β 1 -1 + p 2 Γ(α 2 + β 2 ) Γ(α 2 )Γ(β 2 ) t α 2 -1 (1 -t) β 2 -1 .
In this section, we used the parameters α 1 = 2, β 1 = 5, α 2 = 10, β 2 = 2, p 1 = 2 3 . Graphically, we observe that N (f 1 ) = 2, which is theoretically confirmed in Section 3.2.4.

The Gaussian mixture we chose is the asymmetric claw density introduced by Marron and Wand (1992) :

X ∼    N (0, 1) with probability 1 2 , N l + 1 2 , 2 -l 10 2 with probability 2 1-l 31 , for l ∈ {-2, -1, 0, 1, 2} . (3.13)
Despite the fact that this mixture has 6 components, the underlying density has 5 modes (see [START_REF] Minnotte | The Bumpy Road to the Mode Forest[END_REF]). Its expression is:

f 2 (t) = 1 2 √ 2π e -1 2 t 2 + 2 l=-2 20 31 √ 2π e -50 t-l-1 2 2 -l 2 .
The last mixture studied in this section is made of two Student's t-distributions and satisfies:

X ∼ T (2, -3) with probability 1 2 , T (2, 3) with probability 1 2 , (3.14)
where X ∼ T (ν, ρ) means that X is equal to ρ plus a random variable generated from a t-distribution with ν degrees of freedom. It can be proven that the density f 3 which is associated with the model (3.14) has two modes. In addition, we have

f 3 (t) = 1 4 √ 2 1 + (t + 3) 2 2 -3/2 + 1 4 √ 2 1 + (t -3) 2 2 -3/2
.

Because most of theoretical results we present in this paper concern the L 1 distance between an estimator f and the true density f , it makes sense to use the following criterion, often called integrated absolute error (IAE), defined as:

IAE = f -f L 1 = R f (t) -f (t) dt.
Given different datasets, it can be difficult to compare their respective estimated densities and to draw conclusions from these functions. Locations of modes and antimodes can thus be studied to build real indicators which are easiest to handle as explained in Section 3.2.5. Estimating such positions is straightforward when f and its estimate have the same number of modes. When it is not the case, we would like to determine if these positions can still be estimated. To do so, we estimate the location of the lowest local minimum z = arg min

z i ∈{z 2j } j∈{1,...,N (f )-1} f (z i ),
with z 2j defined in (H3). The estimator ẑ of z is chosen to satisfy ẑ = arg min

z i ∈ Ẑ f (z i ),
where Ẑ is made of the points z such that,

∃ε > 0, ∀ε < ε, f (z -ε) > f (z) and f (z + ε) > f (z).
In , we will study values of IAE as well as estimates ẑ for the three models which are presented in this section and for several density estimators.

Simulation results for the beta mixture

In Figure 3.6, we first draw an example of density estimation with the four considered estimators, using a sample of size n = 1600 generated from model (3.12). Apart from fK,h crit with the uniform kernel, estimates of the density seem to be close to f , considering the shape. We observe a lot of peaks for fP . This feature is directly related to the estimation method as it has been already noticed by [START_REF] Müller | Excess mass estimates and tests for multimodality[END_REF]. Polonik (1995a) wrote that N ( fP ) can be different from N ( f ). In Figure 3.6(b), the estimation of z related to fP is close to 0 and far from z ≈ 0.6219. For a given estimation f of f ,

finding ẑ in ] -∞, z 1 [∪]z 2N (f )-1 , ∞[ only happens if
there is a mode of f in this interval that is sufficiently far from z 1 and from z 2N (f )-1 . This event is especially likely to arise when N ( f ) > N (f ) but it can also happen when N ( f ) = N (f ) if the modes of f are not clearly separated. In addition, when K is the Gaussian kernel, estimations of z for fK,h crit (Figure 3.6(a)) and for fK,h SJ (Figure 3.6(b)) are close to z.

Then, we generate 100 replicates from model (3.12) for various sample sizes n ∈ {100 × 2 i } i∈{0,...,9} . For each sample and each density estimation procedure, IAE is computed. The corresponding values are represented in Figure 3.7 with boxplots. Not surprisingly, for fK,h crit with the uniform kernel, we observe in Figure 3.7(a) that IAE increases with n, which is compatible with Theorem 7. Performances of fK,h crit for the Gaussian kernel, shown in Figure 3.7(b), are better. In this case, boxplots exhibit the L 1 convergence of Theorem 6. They reach a similar precision to those obtained for fK,h SJ which are drawn in Figure 3.7(c). Polonik's method (Figure 3.7(d)) needs extensive computational time. That is why we do not draw boxplots for the greatest sample sizes, but we still observe convergence of this procedure despite of the many peaks of the estimates. Values of IAE appear to be slightly greater for this estimator than those for fK,h SJ and fK,h crit with the Gaussian kernel.

For each replicate we made, we also compute estimations of z. In Figure 3.8 we draw various values of ẑ for the four previously considered estimators. For fK,h crit with the uniform kernel, in Figure 3.8(a), ẑ values move away from the position of the local minimum of f 1 (which is equal to 0.6219) when sample size increases. This is not surprising because of the poor quality of this estimator which, according to Figure 3.6, provides two modes close to each other, but far from the true ones. In Figure 3.8(b), we notice a convergence of the boxplots toward the aimed location for fK,h crit with the Gaussian kernel. For fK,h SJ with the Gaussian kernel, in Figure 3.8(c), a similar convergence can be observed. For fP , in Figure 3.8(d), values of ẑ seem to tend toward 0 or 1 when sample size increases. A possible explanation of this phenomenon is that a spurious mode close to 0 or 1 is sometimes created by the procedure. Because this mode is located in an interval of t where f 1 (t) is small, the local minimum near the mode is the minimum over all local minima. This occurs in Figure 3.6(b), for example. Vertical lines: positions of the minimum of local minima of each plotted density. In Figure 3.6(b), vertical axis is broken between 2.2 and 2.8, between 3.5 and 5 and between 5.5 and 12.5 to be able to clearly see the shape of all curves.

Simulation results for the asymmetric claw density

The simulation results of Section 3.2.3.2 lead us to focus only on the estimators fK,h SJ and fK,h crit with the Gaussian kernel because both exhibit IAE convergence toward 0 and convergence of ẑ toward z. To study them in more depth, we use the model (3.13).

Moreover, this selection is also made because of the costs in computational time of the different methods. For example, for a given sample of size n = 1600, we measured computational time of the methods we consider in this study, with our Intel Core 2 Quad Q9505 processor. To obtain h crit with a Gaussian kernel we need about 1.6 seconds using the density() R function. Finding h SJ requires 0.004 seconds, with the KernSmooth R package while our R implementation of fP needs 425 seconds to be computed. Our R algorithm finds h crit in 4 seconds for the uniform kernel.

In Figure 3.9(a), we draw another example of density estimation using the estimators fK,h SJ and fK,h crit , and a sample of size n = 1600 generated from model (3.13). We observe very similar results for both procedures even if fK,h crit seems to produce an estimate slightly more precise than fK,h SJ for the estimation of z. We also provide estimates fK,h crit,2 and fK,h crit,7 for this sample, in Figure 3.9(b). On one hand, choosing k = 7 > N (f ) does not change the quality of the estimation, comparing to fK,h crit . On the other hand, setting k = 2 < N (f ) leads to an estimator which smooths the modes of f 2 and thus may produce poor values of IAE. Because the only antimode of fK,h crit,2 does not have any reason to be the lowest one, estimating z from this density can also provide inaccurate results, as in Figure 3.9(b). This explains the constraint k ≥ N (f ) in Theorem 6.

Boxplots presenting IAE values for these methods are drawn in Figure 3.10. The dispersion of the IAE values does not seem to decrease with n in Figure 3.10(a). This result could come from the fact that for the asymmetric claw density, (H2) does not hold and Theorem 6 cannot be applied, but IAE values of Figure 3.10(a) globally decrease. Thus, another practical explanation is that the bandwidth h crit cannot adapt to the various sharpness of the modes of f 2 , while fK,h SJ may compensate this phenomenon by creating a new mode such that N ( fK,h SJ ) > N (f 2 ). Indeed, in Figure 3.10(b), IAE values converge toward 0 like those in Figure 3.7(c).

We notice that for the asymmetric claw density, estimates of z can be located at the extrema of the sample, for both fK,h crit and fK,h SJ , in Figure 3.11. This previously occurs for fP and model (3.12) in Figure 3.8(d). For fK,h crit , in Figure 3.11(a), samples that produce this kind of estimate are sufficiently rare for the estimates to be considered as outliers when building boxplots. In Figure 3.11(b), for fK,h SJ and when n ≥ 6400, every ẑ is in the tails of the estimated distribution. Thus, for the estimation of z, h crit seems to perform better than h SJ despite of the fact that (H2) does not hold for f 2 .

Remark 5. [START_REF] Sheather | A reliable data-based bandwidth selection method for kernel density estimation[END_REF] proposed two versions of their bandwidth h SJ respectively called "direct plug-in" and "solve the equation". In this article, only the "direct plug-in" approach is considered. The "solve the equation" version leads to smaller values of IAE in simulation than the "direct plug-in" but needs more computational time. Results still are of the same order of magnitude. Conclusions about the estimation of z are valid for both of them.

Simulation results for the mixture of t-distributions

Studying the model (3.14), we estimate f 3 with fK,h SJ and fK,h crit from a sample of size n = 1600, and with the Gaussian kernel K. In Figure 3.12, we draw the resulting estimates. The bandwidth h crit seems too large since fK,h crit is smoother than f 3 . Because f 3 has two modes, we have N ( fK,h crit ) = 2 but one mode of fK,h crit is close to the modes of f 3 (Figure 3.12(a)) and the other one is in its tails (Figure 3.12(b)). The estimate fK,h SJ appears to perform better than fK,h crit in this example but it also produces a lot of antimodes in the tails. The lowest local minimum of fK,h SJ is thus far from 0 as presented in Figure 3.12(b). In this paragraph, K is the Gaussian kernel. The tendency in Figure 3.12 is confirmed in Figure 3.13 where values of the integrated absolute error between f 3 and fK,h crit (resp. fK,h SJ ), are displayed. When the sample size increases, fK,h crit exhibited poor results. Note that f 3 has heavy tails and then when n is sufficiently large, there is a great probability that a few points of the studied sample are really far from 0. If this happens, the estimator fK,h crit produces a mode near these points and oversmooths the density estimate. Observing Figure 3.13(b), fK,h SJ seems to perform pretty well even when the density to estimate has heavy tails. When n is large, the bandwidth h SJ becomes however really little, implying many spikes and almost constant parts in fK,h SJ . These features lead us to rely on a Monte-Carlo algorithm to integrate since it performs better in such cases than the QUADPACK library used until now. The poor results of fK,h crit encourage us to consider fP again because it does not need that the support of f is bounded. In Figure 3.13(c), we observe that this estimator produces small values of the IAE, even when the tails of f are heavy.

With the Gaussian kernel K, fK,h SJ (t) is often numerically equal to 0 where t is one of the numerous local minimum of fK,h SJ . For example, for 100 samples of size n = 128000, we found that the mean number of antimodes where fK,h SJ (t) = 0 is equal to 813.56. It is thus impossible to find a relevant estimate of z using h SJ . It is also difficult to estimate z from fK,h crit where K is the Gaussian kernel as one can see in Figure 3.14(a). This graphic shows the increasing dispersion of the position of the unique local minimum of fK,h crit when n increases, which results from the presence of modes far from 0. The bandwidth h crit does not lead then to a satisfying estimator of z in this case. Polonik's estimator also exhibits an increasing dispersion in the corresponding values of ẑ as displayed in Figure 3.14(b). Although this increase is less important than the one related to fK,h crit , fP does not provide accurate estimates of z and should not be used for this purpose. While for model (3.13), the assumption (H2) did not appear to be important, we thus see why it is for this mixture of t-distributions. This is consistent with Hall and York (2001, Section 2.3) in which it is stated that h crit,1 does not satisfies (F4) when f is a t-distribution. The authors also indicate that, when f is a normal distribution, if h crit,1 satisfies (F4) then its rate of convergence toward 0 is slow. This could explain Figure 3.10(a). To avoid difficulties when studying densities with unbounded support, [START_REF] Hall | On the calibration of Silverman's test for multimodality[END_REF] propose to replace the function N in (3.9) with a function N I that returns the number of modes of a density which are observed in a given bounded interval I. This interval must be chosen by the user.

Assuming the number of modes of a mixture density

Estimation of a density of a mixture model with a known maximum number of components is the main task that fK,h crit,k can realize. Indeed, if each component of the mixture model is an unimodal density, there are various cases where the number of modes of the density of the mixture is at most equal to the number of components. However, it is not true for every mixture model. For instance, densities made from components f µ,θ that satisfy

f µ,θ (x) = C θ (1 -(x -µ) 2 ) θ 1 [-µ,µ] (x)
, where µ and θ should be chosen for each component, can have a number of modes greater than its number of components (see Remark 3 and Hall et al. (2004)).

Notice that a function which is convex on an open interval does not have any mode on this interval and that the density of a mixture of densities which are convex on an interval is convex on this interval too. Thus, let f be a density of a mixture model with m components. For i ∈ {1, . . . , m}, let g i be the density of the component i. If there exists

µ i ∈ R such that g i is convex on ] -∞, µ i [ and on ]µ i , ∞[, then N (f ) ≤ m and ∀j ∈ {1, . . . , N (f )} , z 2j-1 ∈ ∪ m i=1 {µ i }.
For example one can take for g i the standard two-sided power distribution of van [START_REF] Van Dorp | A Novel Extension of the Triangular Distribution and Its Parameter Estimation[END_REF]:

g i (t) =    γ i t θ i γ i -1 for t ∈ [0, θ i [, γ i 1-t 1-θ i γ i -1 for t ∈ [θ i , 1],
with θ i ∈ [0, 1] and γ i > 2. If we allow ourselves not to respect (H2), the Laplace distribution is another valid choice for g i .

Mixture model densities with N (f ) ≤ m are not restricted to those that satisfy the previous condition of convexity. Considering the beta mixture model of Section 3.2.3 with α 1 = 2 and β 2 = 2, the density f 4 of (3.12) can be written

f 4 (t) = p 1 g 1 (t) + (1 -p 1 )g 2 (t), with g 1 (t) = Γ(α 1 +β 1 ) Γ(α 1 )Γ(β 1 ) t(1 -t) β 1 -1 and g 2 (t) = Γ(α 1 +β 1 ) Γ(α 1 )Γ(β 1 ) t α 2 -1 (1 -t).
When β 1 > 2, for q ∈ {1, 2, 3}, the q th derivative of g 1 satisfies:

g (q)
1 ≥ 0 ⇔ (-1) q-1 x ≤ (-1) q-1 q β 1 , and we also have for q ∈ {1, 2, 3}:

g (q) 2 ≥ 0 ⇔ x ≤ 1 - q α 2 .
These relations hold for strict inequalities and imply that f 4 has no mode in [0,

1 β 1 [ and in ]1 -1 α 2 , 1] because g (1)
1 and g

(1)

2 have the same sign on these intervals on which this sign is constant. Assume now that 3

β 1 < 1 -2 α 2 and that 2 β 1 < 1 -3 α 2
, which is true for the values of β 1 and α 2 chosen in Section 3.2.3 for f 1 . Then, f 4 has at most one mode on

1 β 1 , min 3 β 1 , 1 -3 α 2
. This can be proven by assuming that f 4 has two modes located in z 1 and z 3 and an antimode located in z 2 with 1

β 1 ≤ z 1 < z 2 < z 3 < min 3 β 1 , 1 -3 α 2 . This implies that f (2) 4 (z 2 ) = p 1 g (2) 1 (z 2 ) + (1 -p 1 )g (2)
2 (z 2 ) ≥ 0, and because g

2 (t) > 0 for t ∈ z 2 , min 3 β 1 , 1 -3 α 2 , p 1 g (2) 1 (t) + (1 -p 1 )g (2) 2 (t) > 0 which negates the fact that f (2) 4 (z 3 ) = p 1 g (2) 1 (z 3 ) + (1 -p 1 )g (2)
2 (z 3 ) ≤ 0. A demonstration of the same type leads to the property that f 4 has at most one mode on max 3

β 1 , 1 -3 α 2 , 1 -1 α 2 .
Thus, in order to demonstrate that N (f 4 ) ≤ 2, we have to show that f 4 has no mode on ∆, with:

∆ = max 3 β 1 , 1 - 3 α 2 , min 3 β 1 , 1 - 3 α 2 .
Because of assumptions we made on

β 1 and α 2 , ∆ is included in 2 β 1 , 1 -2 α 2
on which both g

(2)

1 and g

(2)

2 are positive. This implies that f 4 is convex on this interval and then has no mode on it. Thus f 4 has no mode on ∆ and N (f 4 ) ≤ 2. Although we anticipate that generalizing this result to a wider set of mixtures is feasible, the demonstration would probably be tedious.

Oyster opening amplitudes modeled with a bimodal density

In this section we describe a real data application. We apply the estimator fK,h crit,k to opening amplitudes of oysters. These animals are studied by a laboratory called Environnements et Paléoenvironnements Océaniques et Continentaux (http://molluscan-eye. epoc.u-bordeaux1.fr, EPOC) in order to derive water quality indicators. Their approach is based on the assumption that a water of poor quality leads to perturbations in the oysters' behavior. It consists in measurements of the distance between both parts of the shell of the oysters with a frequency of 0.625 Hz. The procedure carried out to obtain the dataset is non invasive. It relies on electrodes stuck on the shell of the oysters and on the GSM/GPRS service to transfer the data. The animals studied in this article are immersed in the Bay of Arcachon, in France. We aim at estimating the density f of the distances of the parts of the shell of these animals. During a day, following the tide, an oyster is either open or closed (see for instance [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF]). For each of this state, the density of the opening amplitudes is assumed to be unimodal. If we also assume that these densities behave similarly to the two-sided power distribution or to the beta mixture detailed in Section 3.2.4, we have N (f ) ≤ 2. The density f has obviously a bounded support since the measured distance is at least 0 and since the opening of an oyster can not be infinitely large.

Because of these assumptions, we estimate f with fK,h crit,2 with the Gaussian kernel and two samples of respective size n = 53932 and n = 53728, which respectively leads to Figure 3.15(a) and Figure 3.15(b). In Figure 3.15(a), the data come from an oyster that does not exhibit any feature of sickness. In Figure 3.15(b), the oyster analyzed died one week after these measures were recorded. Generally, when they are in the open state, dying oysters produce opening amplitudes with a wide variability which implies a larger mode for the corresponding density. Thus, the local minimum z 2 of f located between its two modes comes close to the location z 1 of the mode related to the close state, when oysters' health becomes poor. This feature can also be observed for fK,h crit,2 in Figure 3.15. Indeed, let ẑ1 and ẑ3 be the local maxima of fK,h crit,2 and let ẑ2 be its local minimum such that ẑ1 < ẑ2 < ẑ3 . Then, ẑ2 -ẑ 1 ẑ3 -ẑ 1 is greater in Figure 3.15(a) than in Figure 3.15(b) with respective values 0.3073 and 0.1566. This observation could lead to a detection of oysters in poor health.

Such an indicator can be built with other estimators of f . The point z 2 could be thus estimated by ẑ. Both local maxima for t ≤ ẑ and for t > ẑ could be chosen as estimators of the modes. The simulation study of Section 3.2.3 shows however that ẑ provides good results when it comes from fK,h crit , for densities f with bounded support. This explains the choice of the estimator of f in this section.

Maximum amplitude of the openings of the animal is also linked to its health and we can see that it is approximately 4 times greater in Figure 3.15(a) than in Figure 3.15(b). However, this quantity may also vary with the size of the animal, and with the position of the electrodes on it. That is why we prefer not to rely on it.

Concluding remarks

The estimator fK,h crit,k , when K is a Gaussian kernel, is able to estimate a density f that has a known number of modes because when sample size is large enough, fK,h crit,k is close to f from both L 1 and L ∞ points of view. This is not the case when K is the uniform kernel in our simulation study and theoretical reasons were explored (Theorem 7).

When K is the Gaussian kernel and f has a bounded support, the fact that we require that N ( fK,h crit,k ) = k does not imply a great loss of convergence rate of the IAE for fK,h crit compared with results of fK,h SJ , as shown in our simulation study. In addition, fK,h crit exhibits slightly better results than fP , in this case. When f has heavy tails, fP becomes however more accurate than fK,h crit from the point of view of the IAE. To handle this case, the definition of fK,h crit can be modified such that spurious modes in its tails are ignored (see for instance [START_REF] Hall | On the calibration of Silverman's test for multimodality[END_REF]).

An interesting feature of fK,h crit,k is that it has a deterministic number of modes and can have as many modes as f . In real data analysis, this allows to seek the positions of the various modes and antimodes of the estimate of a density with bounded support, and to derive indicators from them in order to compare densities (see Section 3.2.5). This goal could also be achieved when N (f ) = 2 by fP and fSJ by choosing ẑ as an estimate of the location of the antimode but results for this estimator are better when it is computed from fK,h crit,k , with the Gaussian kernel K. From a practical point of view, note that the bandwidth h crit required more time to compute than h SJ but less than fP .

An important assumption in order for fK,h crit,k to converge toward f is that k ≥ N (f ). In practice, for a mixture model with m unimodal components, this assumption can be replaced by k ≥ m, provided that m ≥ N (f ). Note however that for some mixture models, we have m < N (f ). The estimator fK,h crit,k should thus be used with care in this context.

To conclude, asymptotic properties of fK,h crit with the Gaussian kernel, together with its behavior in simulation, and its deterministic number of modes allow this estimator to be applied to real datasets that are assumed to come from mixture model densities. For that matter, an implementation of this work is included in numerical procedures daily performed on the environmental data of EPOC laboratory.

Microclosings, noise and wavelets

Wavelets are common tools to analyse a signal with various resolutions. They share their theoretical background with Fourier series but they may be more appropriate than these series to describe local behaviours. They are involved in the JPEG-2000 standard for image compressing. Such problems are described in Mallat (2008, Chapter 10). We are here interested in wavelets in order to remove a possible noise from a valvometric signal and to automatically detect spikes in this signal. We focus on Haar wavelets which are piecewise constant because we believe that valvometric signals, without noise, vary slowly in the intervals where there is no spike. We introduce these wavelets in Section 3.3.1. More general kinds of wavelets can be found in [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]. In Section 3.3.2, we explain how to obtain a noise-free signal thanks to wavelets and using especially the works of [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF] and [START_REF] Coifman | Wavelets and Statistics[END_REF]. Among important articles about this subject there are [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF][START_REF] Polonik | Density estimation under qualitative assumptions in higher dimensions[END_REF][START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF], 1999) and Hall and Patil (1996a,b). The locations where a valvometric signal exhibits high variations are likely to correspond to microclosings or transitions between the open and the closed state. In Section 3.3.3, we propose a method which is linked with wavelets in order to find these locations. Walnut (2004, Definition 5.11), describes the Haar wavelets as functions h j,k defined for (j, k) ∈ Z 2 and x ∈ R by

The Haar wavelets

h j,k (x) := 2 j/2 I [ k 2 j , k+1/2 2 j [ (x) -I [ k+1/2 2 j , k+1 2 j [ (x) ,
where for a given subset A ⊂ R, the indicator function I A is equal to 1 when x ∈ A and 0 when x / ∈ A. Let also define for all (j, k) ∈ Z 2 the Haar scaling functions p j,k which satisfies

p j,k (x) := 2 j/2 I [ k 2 j , k+1
2 j [ (x), for every x ∈ R. In this subsection, we only consider functions in L 2 ([0, 1]), which leads us to study only h j,k and p j,k for j ∈ N and k ∈ D j := {0, . . . , 2 j -1}. These functions are interesting because for a given j 0 ∈ N, the indexed family U j 0 := {p j 0 ,k 1 , h j,k 2 } j≥j 0 ,k 1 ∈D j 0 ,k 2 ∈D j is a complete orthonormal system on [0, 1] (Walnut, 2004, Theorem 5.24). This means that U j 0 is an orthonormal basis of the functional space L 2 ([0, 1]). We thus have ∀j 0 ∈ N:

• ∀G ∈ L 2 ([0, 1]), ∀x ∈ [0, 1], G(x) = h∈U j 0 G, h h(x), (3.15)
where for a given (G, h)

∈ L 2 ([0, 1]) × L 2 ([0, 1]), G, h := R G(x)h(x)dx, • ∀(h 1 , h 2 ) ∈ U 2 j 0 , such that h 1 = h 2 , h 1 , h 2 = 0, • ∀h ∈ U j 0 , h, h = h 2 2 = 1.
Equation (3.15) implies that G is fully characterized by U j 0 and by the set { G, h } h∈U j 0 . Assume now that we only have y k = G k 2 J for all k = D J and a given J ∈ N * , at our disposal. Let G J be the piecewise constant approximation of G defined for all x ∈ [0, 1] by:

G J (x) := k∈D J y k I [ k 2 J , k+1 2 J [ (x). If G is continuous, we have G -G J ∞ → 0 and then G -G J 2 → 0, when n → ∞.
Following Walnut (2004, Chapter 6), for all h ∈ {h j,k , p j,k } j∈N,k∈D j ,

G, h -G J , h = G -G J , h ≤ G -G J 2 h 2 = G -G J 2 ,
and then

G J , h → G, h when n → ∞. It thus makes sense to approximate G, h by G J , h for all h ∈ {h j,k , p j,k } j∈N,k∈D j .
Computing these values of G J , h is easy. We indeed have for all k ∈ D J , G J , p J,k = 2 -J/2 y k . We also have for all j > J, for all k 1 ∈ D J and for all k 2 ∈ {2 j-J k 1 , . . .

2 j-J (k 1 + 1) -1}, G J , p j,k 2 = 2 -j/2 y k 1 . In addition, since for all x ∈ [0, 1], G J (x) = k∈D J G J , p J,k p J,k (x),
and because U J is a complete orthonormal system on [0, 1], we have G J , h j,k = 0, for all j ≥ J, k ∈ D j . Only G J , p j,k and G J , h j,k for j ∈ {0, . . . , J -1} and k ∈ D j remain to be determined. Note that for all j ∈ N * and for all k ∈ D j , we have:

G J , p j-1,k = 1 √ 2 G J , p j,2k + 1 √ 2 G J , p j,2k+1 , (3.16)
and

G J , h j-1,k = 1 √ 2 G J , p j,2k - 1 √ 2 G J , p j,2k+1 , (3.17)
as mentioned in Walnut (2004, p142). This allows us to compute recursively G J , p j,k and G J , h j,k for j ∈ {0, . . . J -1} and k ∈ D j , using { G J , p J,k } k∈D J . Thanks to G J , h for h ∈ {h j,k , p j,k } j∈N,k∈D j , we are now able to approximate the decomposition of G into the basis U j 0 , for any j 0 ∈ N. In addition, the quality of this approximation is linked to the closeness between G and G by the Plancherel formula (see for instance Mallat, 2008, Appendix A.3):

G -G J 2 2 = h∈U j 0 G, h -G J , h 2 . (3.18)
Note however that it is not really informative to use the bases U j 0 for j 0 ≥ J.

Denoising with wavelets

Let us now consider the following model for all k ∈ D J :

ỹk = y k + z k , (3.19)
where ỹ0 , . . . , ỹ2 J -1 are available measures, z 0 , . . . , z 2 J -1 are independently and identically distributed from the probability distribution N (0, σ 2 ) and σ 2 is an unknown positive real. For a given j 0 ∈ N, we are interested in estimating G, h for all h ∈ U j 0 . Although wavelet coefficients related to {z k } k∈D J are often considered through a Gaussian white noise we study here the following functions Z J , for the sake of simplicity. They are defined for all J ∈ N * and for all x ∈ [0, 1] by

Z J (x) := k∈D J z k I [ k 2 J , k+1 2 J [ (x).
This definition and equations (3.16) and (3.17) imply that for all j ∈ N, for all k ∈ D j ,

Z J , p j,k ∼ N (0, 2 -J σ 2 ),
and that for all j ∈ {0, . . . , J -1}, for all k ∈ D j ,

Z J , h j,k ∼ N (0, 2 -J σ 2 ).
Using wavelets on G J + Z J , it is thus possible to find an estimate G J of the function G from model (3.19). One may also control G J -G 2 2 thanks to equation (3.18). As explained in [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF], it is expected that for a given j 0 ∈ {0 . . . , J}, only a few coefficients among { G J , h } h∈U j 0 are not null. For any

h 1 ∈ U j 0 such that G J , h 1 = 0 and any h 2 ∈ U j 0 such that G J , h 2 = 0, if σ 2 is small enough, we have with a great probability that | G J + Z J , h 1 | is smaller than | G J + Z J , h 2 |.
It then makes sense for all h ∈ U j 0 to find a threshold δ h > 0 and to choose

G J , h = HT( G J + Z J , h ) := 0 if | G J + Z J , h | < δ h , G J + Z J , h otherwise.
This approach is called hard thresholding. Soft thresholding is another popular method to estimate G, h , for all h ∈ U j 0 , it is defined by:

G J , h = ST( G J + Z J , h ) :=      0 if | G J + Z J , h | < δ h , G J + Z J , h -δ h if G J + Z J , h > δ h , G J + Z J , h + δ h if G J + Z J , h < -δ h .
Both hard and soft thresholding provide then an estimate G J of G, which satisfies for all

x ∈ [0, 1] G J (x) = h∈U j 0 G J , h h(x). (3.20)
To use model (3.19) for valvometric data we have to decide which thresholding method we want to apply, choose j 0 ∈ N to define a basis U j 0 , and set the values of δ h for each h ∈ U j 0 . Let n := 2 J . Various choices are available and have good asymptotic properties such as the SureShrink estimator (Donoho and Johnstone, 1995), the optimal thresholds from a MISE point of view (Hall and Patil, 1996a) or the simple so-called universal threshold defined for all h ∈ U j 0 by δ h = σ 2 log(n)/

√ n, for a given estimator σ of σ [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF] studied various wavelet thresholding estimators in simulation to obtain insights about how they behave when n is finite. They especially consider the case where G is the Bumps function from [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], which is a sum of functions t → (1 -|t|) -4 with various translation and scale parameters. This Bumps function seems to us to be the most interesting one among the functions [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF] tested, because it looks like the valvometric signals we want to denoise. This is mostly due to the spikes this function exhibits. For this function, the best estimator among those considered by [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF] appears to be the translation invariant hard thresholding estimator from [START_REF] Coifman | Wavelets and Statistics[END_REF]. We will detail their approach in Section 3.3.3. Note that we set j 0 = ⌊log(log(n))/ log(2)⌋ + 1 and δ h = σ 2 log(n log(n)/ log(2))/ √ n, where ⌊t⌋ is the integer part of t, for all t ∈ R, like [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: a comparative study[END_REF] did. We discuss a way to find σ in Remark 6.

Finding microclosings with wavelets

Assume that we measure a valvometric signal at each time k 2 J for k ∈ D J . This can be done easily by choosing suitable starting time and unit of time. We can then use the model (3.19) for these data. Microclosings arise for values of t where |G ′ (t)| is great. We also consider that changes between states open and close happen at such instants (see Section 3.1). Provided a way to detect these state changes, the only instants left when |G ′ | is large correspond to the microclosings. We are then interested in studying G ′ . [START_REF] Sow | Water quality assessment by means of HFNI valvometry and high-frequency data modeling[END_REF] estimated G ′ with a kernel-based estimator and did not have to assume that z k comes from a Gaussian distribution, for k ∈ D J . Here, we would like to approximate values of

G ′ k 2 J by 2 J (y k+1 -y k ) for k ∈ D J \{2 J -1} or by 2 J (y k -y k-1 ) for k ∈ D J \{0}. For all even k in D J , we have 2 J (y k+1 -y k ) = 2 3J/2 ( G J , p J,k+1 -G J , p J,k ) = -2 (3J+1)/2 G J , h J-1,k/2 ,
because of equation (3.17). The same equation implies that for all odd k in D J

2 J (y k -y k-1 ) = 2 3J/2 ( G J , p J,k -G J , p J,k-1 ) = -2 (3J+1)/2 G J , h J-1,(k-1)/2 .
Using hard thresholding, and δ h = σ 2 log(n log(n)/ log(2))/ √ n, this leads us, for all

k ∈ D J-1 , to use -2 (3J+1)/2 G J , h J-1,k as an estimator of both G ′ 2k 2 J and G ′ 2k+1 2 J
. Let now assume that we observe {ẙ k := ỹk+1 } k∈D J instead of {ỹ k } k∈D J . This allows us to study G ′ on 2 -J , 1 + 2 -J instead of on [0, 1] by translating the studied functions. We then define GJ and ZJ for all x ∈ [0, 1] by

GJ (x) := k∈D J y k+1 I [ k 2 J , k+1 2 J ] and ZJ (x) := k∈D J z k+1 I [ k 2 J , k+1 2 J ] . For all k ∈ D J-1 , let dk := G J + Z J , h J-1,k = -2 -(J+1)/2 (ỹ 2k+1 -ỹ2k ), and dk := GJ + ZJ , h J-1,k = -2 -(J+1)/2 (ẙ 2k+1 -ẙ2k ).
We then have that -2 (3J+1)/2 HT( dk ) is an estimate of G ′ 2k 2 J and G ′ 2k+1 2 J

, while we also have that -2 (3J+1)/2 HT( dk ) is an estimate of G ′ 2k+1 2 J and G ′ 2k+2 2 J . Both estimates can however take significantly different values. To illustrate this idea, we analyse a valvometric signal of 2 J + 1 = 257 points from an oyster from Eyrac, studied on May 8, 2007 from 00:34:10 to 00:40:59. We display it in Figure 3.16. We can observe a big spike near 0.58 h and another one close to 0.65 h. Smaller ones are located near 0.59 h, 0.62 h and 0.67 h. The corresponding values of {| dk |} k∈D J-1 and {| dk |} k∈D J-1 are plotted in Figure 3.17. They were computed using the wavethresh R package. The hard thresholding sets each value of {| dk |} k∈D J-1 and {| dk |} k∈D J-1 below the dashed line to 0. Both {HT( dk )} k∈D J-1 and {HT( dk )} k∈D J-1 lead to the detection of both main spikes of Figure 3.16. Only {HT( dk )} k∈D J-1 manages to detect the second small spike and only {HT( dk )} k∈D J-1 succeeds in finding the third one. It is then unclear whether or not we should consider that these spikes exist. For denoising purposes this is not an important issue since the inverse wavelet transforms (3.20) respectively computed from {ỹ k } k∈D J and {ẙ k } k∈D J are very similar. We plot them in Figure 3.18.

Notice that for k

∈ D J-1 , 2 J-1 (y 2k+2 -y 2k ) is an approximation of G ′ 2k+1
2 J based on central difference. This encourages us to use -2 (3J-1)/2 (HT( dk ) + HT( dk )) as an estimator of G ′ 2k+1 2 J

. For the signal from Figure 3.16, we obtain the following estimate This vector should be multiplied by 3600 1.6n to be compatible with the scale of Figure 3.16. For each detected spike, the corresponding value -2 (3J-1)/2 (HT( dk ) + HT( dk )) is then coherent with a visual inspection of Figure 3.16.

of G ′ 1 2 J , G ′ 3 2 J . . . , G ′ 2 J -1 2 J ( 13 0
We present now the approach from [START_REF] Coifman | Wavelets and Statistics[END_REF]. For all (k, l) ∈ D 2 J , let ỹk,l be defined by: ỹk,l := ỹ(k+l) mod 2 J , where mod denotes the modulo operator. This means that {ỹ k,l } k∈D J and {ỹ k } k∈D J are linked through a circulant shift. Let also each element of {z k,l } k∈D J satisfy z k,l := z (k+l) mod 2 J , and define G J,l and Z J,l from {ỹ k,l } k∈D J and {z k,l } k∈D J like G J and Z J were defined from {ỹ k } k∈D J and {z k } k∈D J . we will also consider the elements of {d j,k,l } j∈{j 0 ,...,J-1},k∈D j ,l∈D J and {h j,k } j∈{j 0 ,...,J-1},k∈D j where d j,k,l := G J,l + Z J,l , h j,k , and c k,l := G J,l + Z J,l , p j 0 ,k .

One can then use a hard thresholding approach to compute the functions G J,l defined for all x ∈ [0, 1] by

G J,l (x) := k∈D j 0 HT(c k,l )p j 0 ,k (x) + J-1 j=j 0 k∈D j HT(d j,k,l )h j,k (x).
Note that because of the function HT, for each l ∈ D j G J,l changes with j 0 . Coifman and [START_REF] Coifman | Wavelets and Statistics[END_REF] finally propose the translation-invariant estimate G TI,J of G, which satisfies for all x ∈ [0, 1]

G TI,J (x) := 1 2 J l∈D J G J,l x - l 2 J mod 1 .
When j 0 = J -1 and J ≥ 1, we have for all x ∈ [0, 1[

G TI,J (x) = 1 2 G J,0 (x) + G J,1 x - 1 2 J mod 1 .
This is a straightforward consequence of the following property Proposition 8. Let j 0 = J -1 and J ≥ 1. For all (l 1 , l 2 ) ∈ D 2 J such that there exists u ∈ D J-1 which satisfies l 1l 2 = 2u, and for every x ∈ [0, 1[, we have

G J,l 1 x - l 1 2 J mod 1 = G J,l 2 x - l 2 2 J mod 1 .
Proof. This proposition results mainly from the following points:

• for all k ∈ D J-1 \{2 J-1 -1}, for all x ∈ k 2 J-1 , k+1 2 J-1 , h J-1,k (x) = h J-1,k+1 x + 1 2 J-1 and p J-1,k (x) = p J-1,k+1 x + 1 2 J-1 • for all k ∈ D J-1 , for all x ∈ k 2 J-1 , k+1 2 J-1 , for all l ∈ D J , G J,l (x) = HT(c k,l )p J-1,k (x) + HT(d J-1,k,l )h J-1,k (x)
In addition, we have for all k ∈ {0, . . . , 2 J-1 -2} and for all x ∈ k 2

J-1 + 1 2 J , k+1 2 J-1 G J,0 (x) = 2 (J-1)/2 [HT(c k,0 ) -HT(d J-1,k,0 )]
and

G J,1 x - 1 2 J mod 1 = 2 (J-1)/2 [HT(c k,1 ) + HT(d J-1,k,1 )] .
By recalling the definition of dk and dk for all k ∈ {0, . . . , 2 J-1 -2}, we obtain for all

x ∈ k 2 J-1 + 1 2 J , k+1 2 J-1
G TI,J (x) = 2 (J-3)/2 HT(c k,0 ) + HT(c k,1 ) + HT( dk ) -HT( dk ) .

This is coherent with Coifman and Donoho's will to remove spurious spikes from the estimate of G while we use -2 (3J-1)/2 (HT( dk ) + HT( dk )) to detect where each spike is located, whether it is spurious or not.

Remark 6. The translation-invariant thresholding from [START_REF] Coifman | Wavelets and Statistics[END_REF] and our method to find spikes both require an estimate σ of σ because they are based on a threshold δ h which depends on σ. We assume that the valvometric signal is constant when the oyster is closed for a sufficiently large time. It is then possible to estimate σ with the empirical standard deviation σ measured during these periods, but two other problems have to be tackled: for all k ∈ D J , ỹk only takes a finite number of distinct values and besides, the possible electronic noise is more amplified for high values of the signal than for low values. We then consider that for all k ∈ D J , z k follows a binomial distribution with an appropriate scale. For a portion of a valvometric signal recorded when the oyster is closed, let a := min We then set σ = σ c a . We apply our procedure to find spikes in the signal of Figure 2.1. Because we need a sample size n such that it exists a positive integer J satisfying n = 2 J + 1, we focus on the first 32, 769 points of this signal. To estimate σ, we consider the valvometric data between 8.1 and 10.5 hours and the procedure in Remark 6. For 5, 000 points within the studied sample, the results of our detection method is shown in Figure 3.19. If the absolute value of a wavelet coefficient is greater than the threshold δ h , we display it with a cross. We observe that every cross is positioned in the beginning of an important variation of the signal. If a spike is made of a sequence of several great variations, several crosses are then placed on it. Further work should thus be done in order to determine whether or not two consecutive crosses represent the same spike. We could then automatically find the location and the amplitude of each microclosing.

k 1 ∈D J ,k 2 ∈D J {|ỹ k 1 -ỹk 2 | : ỹk 1 = ỹk 2 }.

Implementation

In this section, we describe how the various analyses presented in this chapter were made functional in practice. The main software we use is R. We present it in Section 3.4.1. To handle the valvometric data we studied, we stored them in an SQL database. We explain how to import them into R in Section 3.4.2. Because everyone is not comfortable with the R command-line interface, we also provide a procedure to launch R codes with a graphical user interface in Section 3.4.3.

The R software

Overview. R is both a mathematical software and a programming language. It is a major tool for data analysis, data visualisation and stochastic modelling. It is a free software under General Public License (GPL) and is available for Linux, Mac OS X and Windows. It is based on packages built by the community. A lot of packages and their documentation are stored in the Comprehensive R Archive Network (CRAN, http://cran.r-project.org/). Other package repositories and bundles exist such as Bioconductor (http://www.bioconductor.org/) which contains statistical procedures for biological data analysis. Among commercial alternatives, SAS, Stata and SPSS are commonly chosen by firms.

While the CRAN released informative manuals for both beginners and advanced users, they may find interesting information in books such as [START_REF] Teetor | R Cookbook[END_REF] for quick answers to simple questions and details about tricky pitfalls, or [START_REF] Adler | R in a Nutshell[END_REF] for a more extensive exploration of R. French speakers may also be interested in Lafaye De Micheaux et al.

(2010).

A few notable features. R has many assets we did not need to program our methods. Some of them still deserve to be highlighted. The R language includes a framework to program in an object-oriented way. This framework may be a little difficult to apprehend because, as Adler (2010, Chapter 10) point it out, there exist two different syntaxes to design objects. While new programs should be written with the most recent one (S4), it might be more appropriate to present here the slightly more simple and older one (S3). Using the latter syntax, Algorithm 3.2 creates a function my.function which behaves differently whether its parameter is a real or a character string. Usual functions like print, plot and summary were implemented with S3. If my.function is an S4 generic function, showMethods("my.function") displays R provides tools to visualise multivariate data with static graphics. This is not an easy task when the dimension of the data is greater than 2. The lattice package gathers some of these tools. Algorithm 3.3 underlies the usefulness of this package. This code creates three samples of n = 1000 pairs (X, Y ) which satisfies

   Y = (X -µ) 3 + ε, X ∼ N (µ, 1), ε ∼ N (0, 2.25).
(3.21)

Each sample corresponds to a value of µ in {1, 2, 3}. The random variables X and ε are independent. When the whole samples are pooled, the plot of Y versus X is cluttered and discerning how these variables are linked is difficult (Figure 3.20(a)). The function xyplot can plot separately each sample using the pooled one, if the underlying values of µ are known (Figure 3.20(b)). The function that links Y and X appears then clearly. For each subfigure in Figure 3.20(b), the selected value of µ is given by the grey mark in the strip above it.

How to access MySQL databases with R

Databases are efficient structures to gather various data related to a project. They allow their users to easily retrieve a part of their dataset, using a particular language. SQL (for Structured Query Language) is one of them. The MySQL software is a bundle of tools built for importing data into a database and manipulating them with the SQL language. The mysqlimport program allows us to import valvometric data from text files into a MySQL server every day. This is automatically done with the cron program which can be managed with crontab. We named our database dbValvo and the table that contains the valvometric measures is called tableValvo. The column names of this table can be found using Algorithm 3.8. The French labels lieu, animal, jour, heure and mesure respectively mean "location", "animal", "day", "hour" and "measure".

The RMySQL package Accessing the MySQL command-line interface is not the only way to execute SQL commands. This can also be done with the mysql program and the option -e. In Algorithm 3.9, file.txt receives the result of the SQL command SQL.command. For example, to obtain the valvometric signal of the first oyster of the Eyrac pier, on May 1, 2007, we use Algorithm 3.10.

Analysing data from an SQL database with R and the previous method may be a burden since it requires a text file. Hopefully, the RMySQL package can directly import data from a MySQL server into a data.frame, in a R session. Algorithm 3.11 is a framework to do so. In this algorithm, 0.0.0.0 should be the IP address of the MySQL server, user.name and user.pwd should be valid login password and SQL.command should be an SQL statement.

Installing RMySQL can be quite difficult for Windows users. We provide here a method that worked for a PC with Windows 7 and RMySQL 0.9-3.

• Install R, MySQL, RTools and the DBI package.

• Download RMySQL source files.

• Extract the .tar archive from the .tar.gz archive and then extract the files from the .tar archive into a RMySQL folder.

• Open the Command Prompt. -------+-------------+------+-----+---------+-------+ -------+-------------+------+-----+---------+-------+ 5 rows in set (0.00 sec) Algorithm 3.8: Displaying the column names of a table

> n <-1000 > nMu <-3 > mu <-rep(1:nMu, n) > X <-rnorm(n * nMu, mu) > Y <-(X-
+--------+-------------+------+-----+---------+-------+ | Field | Type | Null | Key | Default | Extra | +-
| lieu | varchar(50) | YES | MUL | Eyrac | | | animal | tinyint(4) | YES | | NULL | | | jour | date | YES | | NULL | | | heure | double | YES | | NULL | | | mesure | double | YES | | NULL | | +-
• Enter: SET CYGWIN=nodosfilewarning.

• Go to the folder where the RMySQL folder is, using the command cd which is available since RTools is installed.

• If R folder is the path to the folder where R is installed, enter "R folder\bin\ R.exe" CMD build RMySQL to build the RMySQL package.

• Click on "Start Menu", "Control Panel", "System", "Advanced System Settings" and "Environment Variables".

• Edit the system variable PATH by adding c:\Rtools\gcc-x.y.z\bin; at the beginning. Use you file explorer to find the appropriate values of x, y and z.

• In a new system variable MYSQL HOME, write the path to the folder where MySQL is installed, like it should be in Linux (for example C:/Program Files/MySQL/MySQL Server 5.5).

• In the folder where MySQL is installed, find the lib folder and create an opt folder in it. Move the file libmysql.lib from lib to opt.

• Back in the Command Prompt, go where the RMySQL folder is.

• Enter: "R folder\bin\R.exe" CMD INSTALL RMySQL x-y-z.tar.gz. Use you file explorer to find the appropriate values of x, y and z.

• Launch R and type library(RMySQL).

Remark 7. We created an index on the table dbValvo.tableValvo in order for the server to return fast answers when an appropriate WHERE clause is provided. The index sort the rows of dbValvo.tableValvo following the column lieu, and then the column animal and finally jour. This means that requesting every measure from a given location and a given animal may be faster than asking for every measure from a given location and a given day. More details about the indexes can be found in the documentation of MySQL.

How to link R with a graphical user interface

Using Algorithm 3.11 implies writing the password in the R console but graphical objects that hide the entered content exist. Some other programs might need a dynamic adjustment to perform well. This can also be realized with a graphical user interface (GUI) as in [START_REF] Mangiarotti | Polynomial search and global modeling: Two algorithms for modeling chaos[END_REF]. Such interfaces can finally be designed to make user-friendly functions.

GTK+ is a set of tools for creating GUIs. It can be linked with R with the RGtk2 package. GUIs can be designed with the Glade software. We recommend to first install Glade, which will also install GTK+ and then to download RGtk2. Glade allows developers to save their GUI in a .xml format. An example is given in Algorithm 3.12. This file describes a GUI with a window called myWindow and a button named myButton.

Algorithm 3.13 shows how to load this GUI with R. The functions gtkWindowGetType() and gtkButtonGetType() initialize the corresponding types of object. The .xml is imported into the variable builder. The function getObject() allows us to reach Figure 3.21: We designed a graphical user interface to analyse valvometric data each graphical object. For a given kind of object, we can then call the functions defined by GTK+, like show(). See the following documentation for more information: https://developer.gnome.org/gtk2/stable/. The function gSignalConnect should be considered to execute R code when a graphical event arise. In Algorithm 3.13, when myButton is clicked, a new window appears.

For the valvometric data of this chapter, we created the GUI shown in Figure 3.21. A user can enter its login and password to connect to the MySQL database. A calendar then becomes visible. After selecting a location, a month and an animal, he can display a signal by clicking on the corresponding day in the calendar. In another tab, he can browse this signal and remove its noise (see Section 3.3). He can also draw estimates f n of the density f introduced in Section 3.1.4.1 for a given oyster and several days. He is then able to produce a 3-dimensional figure (day vs oyster's opening vs value of the estimates), using the rgl package. For an oyster of Eyrac, from May 2007 to August 2007, an example of such a plot is given in Figure 3.22. In this figure, the location of the mode of f n which corresponds to the lowest part of the valvometric signal increases with time. This illustrates the growth of the animal.

Fast computation with R

Because R is an interpreted language, it is computational time consuming. Using matrix calculus instead of loops produces faster programs. For large project development, calling libraries written in other languages like C or Fortran is recommended. The dyn.load R function loads such libraries. We can access functions in a library thanks to .Call. For example, the density function computes kernel density estimators using C code. The QZ algorithm presented in Section 4.4.2.3 can be reached by .Call from R.

In a command-line interpreter, R CMD BATCH can be executed to obtain the results of some R code without launching R. It is then easy to call R from a C program. We can take advantage of computer clusters using the C language. With this method, we then can do it with R too. An alternative is given by the doMPI package.

$ mysql -u user.name -p -h 0.0.0.0 -e "SQL.command" > file.txt Algorithm 3.9: How to execute an SQL command with mysql mysql> SELECT heure, mesure FROM dbValvo.tableValvo -> WHERE lieu='EYRAC' AND animal=1 -> AND jour='2007-05-01';

Algorithm 3.10: An SQL command to extract valvometric data from dbValvo.tableValvo > library(RMySQL) > con <-dbConnect(MySQL(), host="0.0.0.0", + user="user.name", password="user.pwd", dbname="dbEels") > output <-dbSendQuery(con, "SQL.command") > result <-fetch(output, n=-1) > print(result) ... > dbDisconnect(con) Algorithm 3.11: Sending an SQL command from R Figure 3.22: Estimates of the probability density f for an oyster at the Eyrac pier from May 1, 2007 to August 31, 2007; x: days, y: distance between shell parts, z: estimated density <?xml version="1.0"?> <interface> <requires lib="gtk+" version="2.16"/> <object class="GtkWindow" id="myWindow"> <child> <object class="GtkHButtonBox" id="hbuttonbox1"> <property name="visible">True</property> <child> <object class="GtkButton" id="myButton"> <property name="label" translatable="yes"> Create a new window</property> <property name="visible">True</property> <property name="can focus">True</property> <property name="receives default">True</property> </object> <packing> <property name="expand">False</property> <property name="fill">False</property> <property name="position">0</property> </packing> </child> </object> </child> </object> </interface> Algorithm 3.12: File ex.xml which contains characteristics of a graphical interface > library(RGtk2) > gtkWindowGetType() > gtkButtonGetType() > builder <-gtkBuilder() > builder$addFromFile("ex.xml") > builder$connectSignals(NULL) > myButton <-builder$getObject("myButton") > myWindow <-builder$getObject("myWindow") > myWindow$show() > > gSignalConnect(myButton, "clicked", function(btn){ > myWindow2 <-gtkWindow() > }) Algorithm 3.13: A window which can create new windows Chapter 4

Statistical methods to analyse genetic data Among datasets that help to understand how genetics affect biological observations we will mostly consider gene expressions in this chapter. More specifically, we will study which ribonucleic acid (RNA) sequences are present in a given cell. Such a science is called transcriptomics, because deoxyribonucleic acid (DNA) is said to be transcribed into RNA. While DNA contains information to build proteins, only a part of it is considered at any given moment. A gene is a portion of the DNA. For two distinct genes, we obtain two distinct RNA sequences. An expressed gene is then a gene which corresponds to a RNA sequence we observe.

In Section 4.1, we detail two examples in which gene expressions are involved. Focusing on the second one, we describe in Section 4.2 a method to link gene expressions to an explanatory variable, using linear functions and multiple tests. We propose a Multivariate SIR (MSIR) method based on Sliced Inverse Regression (SIR), which can detect non-linear relationships between two sets of variables, in Section 4.3. We apply it to pupil rating and to a remote sensing dataset, although the MSIR method seems to be able to handle transcriptomic datasets. In Section 4.4 we finally tackle the case when the sample size is smaller than the number of explanatory variables by comparing four improvements of SIR built to solve this problem.

Challenges

Finding expression Quantitative Trait Loci (eQTL)

Concerning diseases which directly result from precise alleles of genes, an important goal is to locate the parts of the genome which are implicated. These illnesses are traits and they are called quantitative traits if they can be quantitatively measured. The portions of DNA which affect the values of these traits are then Quantitative Trait Loci (QTL). Other diseases may occur when a set of genes is too much or too little expressed. See [START_REF] Schadt | Genetics of gene expression surveyed in maize, mouse and man[END_REF] and references therein for examples. While knowing which are these genes is informative, finding those which control how the first ones are expressed is also of particular interest. Such genomic regions are named expression Quantitative Trait Loci (eQTL).

In Section 4.4.5, we try to find eQTL for the Hopx gene in rats. According to WikiGenes (http://www.wikigenes.org) and to [START_REF] Yamaguchi | Homeobox gene HOPX is epigenetically silenced in human uterine endometrial cancer and suppresses estrogen-stimulated proliferation of cancer cells by inhibiting serum response factor[END_REF], Hopx is often called HOP, LAGY or NECC1. This gene allows an organism to create the HOP protein which may be needed for its heart to develop (see [START_REF] Shin | Modulation of cardiac growth and development by HOP, an unusual homeodomain protein[END_REF]. Variations in the expression of this gene can be linked with various types of cancer as in [START_REF] Chen | Identification of a novel homeobox-containing gene, lagy, which is downregulated in lung cancer[END_REF]. That is why, finding eQTL for the Hopx gene seems to be a significant step toward the creation of a therapy for these diseases.

Selecting genes to build an RNA microarray

While in the previous subsection, we focused on finding relevant explanatory variables according to the expression of a given gene, we would like now to keep only genes which exhibit variations in their expression because of some explanatory variable. For example, this need can arise when we have at our disposal gene expressions which come from a DNA sequencing procedure. This technique gives information about every expressed gene of an individual. Because DNA sequencing is an expensive technology, we could prefer microarrays, which are cheaper. Such devices can only measure the expression of a limited number of genes. When the DNA sequencing produces data for a lot of different genes, a selection should then be made.

In the IMMORTEEL project which is supported by both the French ANR and the Canadian CRSNG and is coordinated by Magalie Baudrimont and Patrice Couture, the biologists are facing this problem. They would like to characterise how various pollutants affect eels which spend time either in Canadian rivers or in French ones. For this purpose, gene expressions were measured from these animals, using DNA sequencing. Biologists found 17,866 expressed genes. In her work to obtain a PhD, Lucie Baillon, supervised by Fabien Pierron, would like to build a microarray with only 1,000 genes. This selection should be realized by taking some explanatory variables into account.

These variables aim either at determining how grown the animal is, characterising normal variation of its environment or measuring how much the eel is subject to pollutants. We are especially interested in the last group of variables. For each pollutant, we would like to select genes that are differently expressed whether we detect it in the animal or not. Such genes should correspond to a single pollutant. We thus be able to detect its presence using only gene expressions. The variables which are not related to any pollutant are included in the study in order to select some gene expressions that do not vary according to the standard evolution of the eel or of its environment, but only with a single pollutant. Concerning the growth of the animal we have

• the proportion Lip of lipid in muscles,

• the length L of the animal,

• the relative condition factor (Kn) which indicates how much a fish is heavy with respect to its length,

• the hepatosomatic index (HSI ) which is defined as the ratio between the liver weight and the total weight,

• the splenosomatic index (SSI ) which is the spleen weight divided by the total weight.

The following variables were measured in the fishing sites:

• the salinity S of the water in permille,

• the temperature T of the water inC.

The pollutants taken into account are enumerated in the next list. Metals and arsenic are studied through their weight by gram of dry kidney. The other pollutants are measured in weight unit by gram of muscle.

• Cadmium, copper, zinc, silver, arsenic, lead, chromium, nickel and mercury respectively correspond to the variables Cd, Cu, Zn, Ag, As, Pb, Cr, Ni and Hg,

• several polychlorinated biphenyls (PCBs) are merged in the variable PCB,

• a dichlorodiphenyldichloroethane (DDD) produces the variable DDD,

• a dichlorodiphenyldichloroethylene (DDE) and dieldrin form the variable DDE Diel,

• another DDE leads to the variable DDE,

• a dichlorodiphenyltrichloroethane (DDT) provides the variable DDT,

• the variable Ld is built from lindane,

• the variable PBDE stands for a tetrabromodiphenyl ether called PBDE-47,

• the variable HCB represents hexachlorobenzene (HCB),

• the variable Tnc symbolizes trans-nonachlor.

The considered metals and arsenic can be gathered together in the heavy metal group. This family of elements does not have any consensual definition. It generally contains metals and metalloids that are toxic for some living beings. Cadmium, for example, is responsible for renal problems as explained in [START_REF] Bernard | Cadmium & its adverse effects on human health[END_REF] or in [START_REF] Järup | Health effects of cadmium exposure-a review of the literature and a risk estimate[END_REF].

Although too much copper, zinc or chromium can be harmful for humans, they are also essential for their organism to work properly. Ingestion of lead can cause saturnism. Mercury provokes hydrargyria. The toxicity of a compound can be measured using the LD 50 indicator which is the quantity such that half a group of animals will die if they have to handle it. It is measured in weight of the compound per unit of mass of the animal. Hughes (2002, Table 1) gives the LD 50 for several molecules that contain arsenic. For example, arsenic trioxide (As 2 O 3 ) has an LD 50 of 26 mg/kg for rats exposed by oral route. Silver is not very dangerous for human but the possible effect it could have on aquatic animals is still examined [START_REF] Ratte | Bioaccumulation and toxicity of silver compounds: A review[END_REF][START_REF] Call | Silver toxicity to chironomus tentans in two freshwater sediments[END_REF]. Nickel is studied for its possible carcinogen behaviour and for the inflammations of the skin it may induce. Some toxic metals are not easily excreted by living beings and can stay a long time in their organism. If an animal lives in an environment which contains a low concentration of such a metal, this element can thus affect the animal after a sufficient long time. This is called bioaccumulation. This phenomenon is for example studied in [START_REF] Casas | Modélisation de la bioaccumulation de métaux traces (Hg, Cd, Pb, Cu et Zn) chez la moule, Mytilus Galloprovincialis, en milieu méditerranéen[END_REF]

A regression model with factors 4.2.1 Testing hypotheses with the False Discovery Rate

Assume that we want to study the expression of q genes and that we have n measures for each of them. For j ∈ {1, . . . , q}, let y (j) be the random vector of length n which produces the expressions of the j th gene. Let X = (X 1 , . . . , X p ) ′ be the p × n matrix made of n values of each of the p explanatory variables. We suppose that n > p. For all j ∈ {1, . . . , q}, the following linear regression model may be chosen to link y (j) to X,

y (j) = X ′ β (j) + ε (j) , (4.1)
where β (j) is a real vector of length p and ε (j) is a random vector of length n. For the sake of simplicity, we suppose that ε (j) ∼ N n (0, σ 2 j I n ), where I n is the identity matrix and σ j is an unknown positive real.

Let H 0 be a boolean, and let U be a dataset which comes from a model M (H 0 ). Properties about M (1) and M (0) allows researchers to build statistical tests. A statistical test is a procedure that takes U and a parameter α ∈]0, 1] as inputs and produces

• an estimate H 0 of H 0 ,
• a property which relies on α.

This property is an insight about how much H 0 and H 0 are likely to be the same. Usual statistical tests provide the following property:

P( H 0 = 0|H 0 = 1) ≤ α. (4.2)
Note that these tests also satisfies for all α ∈]0, 1],

P( H 0 = 0|H 0 = 0) → 1, (4.3)
when the sample size n of U goes to ∞. That is why, when H 0 = 0 and α is small enough, we can be confident about the fact that H 0 = 0. For j ∈ 1, . . . , q, if β (j) i

is the i th element of the vector β (j) in (4.1), there exist usual statistical tests for the following hypothesis

H 0 = H (j) 0,i := 0 if β (j) i = 0, 1 otherwise,
with U = U (j) := (y (j) , X). Among them, we choose the F-Test which relies on the property that for any j ∈ {1, . . . , q} and for any i ∈ {1, . . . , p}, if H (j) 0,i = 1:

(np) Hy (j) -H -i y (j) 2 y (j) -Hy (j) 2 ∼ F(1, np),

where

H := X ′ (XX ′ ) -1 X, H -i := X ′ -i (X -i X ′ -i ) -1 X -i , X -i := (X 1 , . . . , X i-1 , X i+1 , X p ) ′ ,
F is the Fisher-Snedecor distribution, and for any vector a, a 2 = a ′ a. We call this test W i . More details about it can be found in Bercu and Chafaï (2007, Chapter 7). If X i 2 → ∞ when n → ∞, W i satisfies equation (4.3). For all α ∈]0, 1], the estimate of H (j) 0,i provided by W i (U (j) , α) is named H (j) 0,i . Let U ⋆ := U (1) , . . . , U (q) , and let W ⋆ i be a sequence of tests defined for a given γ ∈]0, 1] by 1) , γ), . . . , W i (U (q) , γ)}. The False Discovery Rate (FDR) is formalised in [START_REF] Benjamini | Controlling the False Discovery Rate: A practical and powerful approach to multiple testing[END_REF]. It is the expected proportion of false positive tests among positive tests and then satisfies

W ⋆ i (U ⋆ , γ) := {W i (U (
F DR(W ⋆ i (U ⋆ , γ)) := E #{ H (j) 0,i = 0 and H (j) 0,i = 1} #{ H (j) 0,i = 0} .
We consider here that 0/0 = 0. For all j ∈ {1, . . . , q}, let α(j) i be the p-value of W i (U (j) , .), that is α(j) i := sup{α ∈]0, 1] :

H (j) 0,i = 1}.
Without any loss of generality, assume that α(1) i ≤ . . . , ≤ α(q) i . From a usual statistical test W i , and a dataset U i := (W i , U ⋆ ), provided that for all (j 1 , j 2 ) ∈ {1, . . . , q -1}×{j 1 +1, q}, α(j 1 ) i and α(j 2 ) i are independent, [START_REF] Benjamini | Controlling the False Discovery Rate: A practical and powerful approach to multiple testing[END_REF] managed to create a sequence of tests W ( U i , γ) defined by

• estimates H (j) 0,i of H 0,i for all j ∈ {1, . . . , q} satisfying

H (j) 0,i := 0 if ∃j ⋆ ≥ j, α(j ⋆ ) i ≤ j ⋆ q γ, 1 otherwise, • the property F DR( W ( U i , γ)) ≤ γ.
Note that for the estimate H (j) 0,i produced by W i (U (j) , α), we have

H (j) 0,i = 0 if α(j) i ≤ α, 1 otherwise.
For all j ∈ {1, . . . , q}, we can also define a p-value γ(j)

i of W i ( U i , .) by γ(j) i := sup{γ ∈]0, 1] : H (j) 0,i = 1}.
If for a given i ∈ {1, . . . , p}, W i satisfies (4.3), we also have for all j ∈ {1, . . . , q}, and all γ ∈]0, 1] P( H

(j) 0,i = 0|H (j) 0,i = 0) → 1. (4.4)
This can be shown by noticing that if α = 1 q γ, ( H

(j) 0,i = 0) ⇒ ( H (j) 0,i = 0) and if α = j ⋆ q γ, ( H (j) 0,i = 0) ⇒ ( H (j) 0,i = 0).

FAMT: a method to take the correlations between the error terms into account

For a given i ∈ {1, . . . , p} and a given j ∈ {1, . . . , q}, α(j) i relies on ε (j) so that for all (j 1 , j 2 ) ∈ {1, . . . , q -1} × {j 1 + 1, q}, α(j 1 ) i and α(j 2 ) i are independent if and only if the covariance matrix of (ε (j 1 )′ , ε (j 2 )′ ) ′ is diagonal, which means that the link that may exist between y (j 1 ) and y (j 2 ) is fully characterized by X ′ β (j 1 ) and X ′ β (j 2 ) . We would like to weaken this assumption. To do so, like [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF] did, we introduce the following model, for j ∈ {1, . . . , q}

y (j) = X ′ β (j) + Z ′ b (j) + ε (j) , (4.5)
where b (j) is a vector of length r and Z is a r × n random matrix. In addition

, if Z = [Z 1 , . . . , Z n ], we assume that    Z 1 . . . Z n    ∼ N nr (0, I nr ).
We also suppose that the covariance matrix of (ε (1)′ , . . . , ε (q)′ ) ′ is diagonal, but this does not imply this time that for all (j 1 , j 2 ) ∈ {1, . . . , q -1} × {j 1 + 1, q}, y (j 1 ) and y (j 2 ) are only linked through X. [START_REF] Friguet | A factor model approach to multiple testing under dependence[END_REF] use an EM algorithm to find estimates Z and b(1) , . . . , b(q) of Z and b (1) , . . . , b (q) . Their procedure includes a criterion to choose a suitable value of r. Note that for all j ∈ {1, . . . , q}, if y (j) satisfies model (4.5), then ỹ(j) = y (j) -Z ′ b (j) satisfies model (4.1). Let ŷ(j) = y (j) -Z ′ b(j) . For a given i ∈ {1, . . . , p} and any γ ∈]0, 1], we can then perform the sequence of tests W ( U i , γ) using U (j) = ( ŷ(j) , X) for all j ∈ {1, . . . , q}.

Gene expressions of eels studied with FAMT

We analyse the dataset presented in Section 4.1.2. Let X 1 be a vector made of n values of 1. When p ≥ 2, using model (4.5), we would like to have

β (j) = (β (j) 1 , i-2 0, . . . , 0, β (j) i , p-i 0, . . . , 0) ′ , (4.6)
where β (j) i = 0, in order to link a single explanatory variable i ∈ {2, . . . , p} to a gene j ∈ {1, . . . , q}. FAMT allows us to test if β (j) i

is not null but we may not be able to confidently assure that for a given i and a given j, β

(j) i = 0. We nevertheless set β(j) i = 0 if γ(j)
i is greater than a chosen value γ ⋆ , because of equation (4.4) for the sequence of tests W ( U i , γ ⋆ ).

For all j ∈ {1, . . . , q}, and a pair (i 1 , i 2 ) ∈ {2, . . . , p-1}×{i 1 +1, . . . , p}, if cor(X i 1 , X i 2 ) 2 is close to one, it is very unlikely to have γ(j) i 1 ≤ γ and γ(j) i 2 ≥ γ ⋆ for values γ and γ ⋆ such that γ ⋆γ is large. We would like then to withdraw the components of X which are too much correlated with other components. To do so, we display in Figure 4.1 the projection of the explanatory variables onto the first two principal components computed from X. The respective proportions of the total empirical variance of X along these axes are equal We have n = 44 eels. Gene expression data and observations of the explanatory variables were stored in a MySQL server. We access it from R like we did in Section 3.4 for the valvometric signals. We then study it with the FAMT model thanks to the FAMT package [START_REF] Causeur | Factor analysis for multiple testing (FAMT): an R package for large-scale significance testing under dependence[END_REF]. We created a GUI to perform easily such analyses. There are three tabs in this interface. The first one (see Figure 4.2) includes a frame to connect to the SQL server and another one to select the explanatory variables to consider. We can launch the modelFAMT function from this tab and obtain γ(j) i for (i, j) ∈ {2, . . . , p} × {1, . . . , q}. These p-values can be saved from the second tab of our interface (see Figure 4.3). Finally, Figure 4.4 shows the last tab of the GUI which lists genes j such that for a variable X i 1 , γ(j) i 1 ≤ γ and γ(j) i 2 ≥ γ ⋆ for all i 2 ∈ {2, . . . , p}\{i 1 }. For example, setting γ = 10 -3 and γ ⋆ = 0, when X i 1 corresponds to Cd leads to select 215 genes. If we take γ ⋆ = 2 × 10 -3 , we only obtain 92 genes which is a more suitable number since we have space for 1, 000 genes in the microarray and we have 10 pollutants left in the set of explanatory variables. Because γ is already small, we prefer to introduce a non-null γ ⋆ rather than choosing a smaller value for γ.

Introduction

In analyzing large datasets, multivariate response regression analysis with a p-dimensional vector of regressors has been extensively studied in the literature. The reduction of the dimension of the regressors' space is a major concern in this framework. When the response variable is univariate, the issue has been addressed by [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] via the notion of EDR (effective dimension reduction) space. The EDR directions (which form a basis of this subspace) are used to project the p-dimensional covariate x on a K-dimensional linear subspace (with K < p) first for displaying and then for studying its relationship with the response variable y. When the dimension of y is one, it is easy to view the link between the projected predictors and the response variable. The notion of EDR space was also clarified by Cook and his collaborators in their numerous papers introducing the notions of central subspace and central mean subspace, see for details Cook (1998) or [START_REF] Cook | Dimension reduction for conditional mean in regression[END_REF]. [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] introduced sliced inverse regression (SIR) which is a well-known method to estimate the EDR space. The link function can be estimated with a smoothing method such as kernel or smoothing splines approaches for instance.

In this paper, a q-dimensional response variable y is considered. Hence, we deal with a high dimensional regression framework which is not a linear one or a prespecified parametric one. The underlying idea of the dimension reduction of the explanatory variable x without loss of information is to identify the linear combinations

β ′ 1 x, . . . , β ′ K x such that y ⊥ x|(β ′ 1 x, . . . , β ′ K x), (4.7)
where ⊥ denotes independence, K(≤ p) is as small as possible and the p-dimensional vectors β k are linearly independent. Let B = [β 1 , . . . , β K ] denote the p × K matrix of the β k 's. Statement (4.7) means that y|x and y|B ′ x share the same distribution for all values of x. A straightforward consequence is that the p-dimensional covariate x can be replaced by the K-dimensional predictor B ′ x without loss of regression information. The goal of dimension reduction is achieved for K < p. As mentioned in [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] or [START_REF] Cook | On the interpretation of regression plots[END_REF], statement (4.7) is equivalent to y ⊥ x|P B x, where P B denotes the projection operator on Span(B) which is the linear subspace of R p spanned by the columns of B. In addition, Span(B) can be viewed as the EDR space. From a regression model point of view, one can mention that the corresponding underlying model is the following semiparametric one

y = f (B ′ x, ε), (4.8)
where f : R K+r → R q is an arbitrary and unknown link function, ε is a r-dimensional random error variable independent of x (with r ≥ 1). [START_REF] Li | Dimension reduction for multivariate response data[END_REF] consider a regression model with an additive error term: y = g(B ′ x) + ε, where g is an unknown link function taking its values in R q .

In the following, we consider a slightly more restrictive regression model defined as follows:

     y (1) = f 1 (B ′ x, ε (1) ) . . . y (q) = f q (B ′ x, ε (q) ) (4.9)
where for j = 1, . . . , q, y (j) (resp. ε (j) ) stands for the jth component of y (resp. of ε) and the link function f j is an unknown real-valued function. In this paper, we propose an approach to estimate the corresponding EDR space. It is based on combining information from the marginal regression of each component y (j) of y.

Remark 8. The information from the marginal regression is sufficient to recover the whole EDR space in model (4.9). However, this is not always the case when working with model (4.8). Let us illustrate this point with the following regression model proposed by [START_REF] Zhu | On dimension reduction in regressions with multivariate responses[END_REF]:

y (1) y (2) ∼ N 2 0 0 , 1 sin(B ′ x) sin(B ′ x) 1 .
In this example, since the information of interest is only in the correlations between the components of the response variables, none of the available marginal regression provides useful information about the EDR space while considering the entire y allows to recover it.

The vectors β k are not individually identifiable neither in model (4.8), nor in the more restrictive one (4.9). Thus, the main objective is to estimate a basis of the K-dimensional EDR space. Many numerical methods have been introduced to achieve this goal. Let us mention three of them which are relatively simple and easy to implement: SIR introduced by [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] for the single index model (K = 1) and [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] for the multiple indices model (K > 1), principal Hessian directions (see for instance [START_REF] Carroll | Measurement error regression with unknown link: dimension reduction and data visualization[END_REF] or Cook (1998)), and sliced average variance estimation (see for details Cook (2000), Prendergast (2007) or [START_REF] Shao | Partial central subspace and sliced average variance estimation[END_REF]). In this paper, we shall only focus on the SIR approach which is based on a property of the first moment of the inverse distribution of x given y.

From a theoretical point of view, the following linearity condition is required:

(LC) the conditional expectation E[b ′ x|B ′ x] is linear in B ′ x for any b ∈ R p .
One can observe that the linearity condition does not involve the response variable y and only concerns the distribution of the covariate x. Let us mention that, when the distribution of x is an elliptically symmetric distribution (such as a p-dimensional normal distribution), this condition is satisfied. [START_REF] Cook | Reweighting to achieve elliptically contoured covariates in regression[END_REF] proposed a method based on the minimum volume ellipsoid to transform and weight the predictors in order to approximate ellipticity. Kuentz and Saracco (2010) recommended to cluster the predictor space so that the linearity condition approximately holds in the different partitions. To conclude this brief discussion on the linearity condition, using a Bayesian argument of [START_REF] Hall | On almost linearity of low dimensional projections from high dimensional data[END_REF], we can infer that (LC) approximately holds for many high-dimensional datasets (that is when p is large).

As previously mentioned, the proposed approach to estimate the EDR space relies on combining estimates from the marginal regressions of model (4.9). Moreover, we can naturally take the information from these regressions into account to detect if a common EDR space really exists for all the components of y as in model (4.9). Otherwise, one shall consider the following more general regression model for multivariate response regression:

                           y (1) = f 1 (B ′ 1 x, ε (1) ), . . . = . . . y (q 1 ) = f q 1 (B ′ 1 x, ε (q 1 ) ), y (q 1 +1) = f q 1 +1 (B ′ 2 x, ε (q 1 +1) ), . . . = . . . y (q 1 +q 2 ) = f q 1 +q 2 (B ′ 2 x, ε (q 1 +q 2 ) ), . . . = . . . y (q) = f q (B ′ L x, ε (q) ), (4.10)
where, for every l = 1, . . . , L, B l is a p×K matrix, K is assumed to be known, Span(B 1 ) = Span(B 2 ) = . . . = Span(B L ) and L l=1 q l = q. This means that writing this model as in (4.9) requires the number of columns of B to be greater than K. When trying to reduce as much as possible the dimension of a model, estimating B 1 , . . . , B L seems thus more appropriate than seeking B. One then needs to cluster the components of y associated with the same EDR space. Therefore, for each identified cluster of components, one can use only these components to estimate the corresponding (common) EDR space. In a more general case, one can also assume that the dimension K is specific for each B l .

The goal of this paper is twofold. First we propose a new multivariate SIR approach for estimating the K-dimensional EDR space which is common to the q components of the multivariate response variable in model (4.9). Then we propose a way to cluster the components of y associated with the same EDR space in model (4.10) in order to apply properly our multivariate SIR on each cluster instead of blindly applying it on all the components of y.

The paper is organized as follows. In Section 4.3.2, we give a brief overview on usual univariate SIR and existing multivariate SIR methods. The population version of the new SIR approach for a multivariate response, named MSIR hereafter, is described in Section 4. 3.3.1. The corresponding sample version is introduced in Section 4.3.3.2 and asymptotic results are provided in Section 4. 3.3.3. A weighted version of MSIR, named wMSIR hereafter, is proposed in Section 4. 3.3.4. Both these methods rely on a tuning parameter H which is the number of slices. The choice of H and of the dimension K is discussed in Section 4. 3.3.5. Practical methods to investigate the possible existence of a common EDR space for y and to detect and identify clusters of components of y are proposed in Section 4. 3.4. Numerical results based on simulations are exhibited in Section 4.3.5 in order to show the good behavior of MSIR and wMSIR approaches and the usefulness of the diagnostic and clustering procedures on the components of y. In Section 4.3.6, two real datasets are considered: the first one concerns hyperspectral remote sensing while the second one is the widely studied Minneapolis elementary schools dataset. Finally, concluding remarks are given in Section 4.3.7.

Brief review of univariate and multivariate SIR approaches

In this section, we consider the regression model (4.9). We first provide an overview of the SIR method when the response y is univariate. Then, some existing SIR methods for a multivariate response are briefly described. The aim of all these approaches is to estimate the EDR space.

Univariate SIR

We focus here on a univariate response (i.e. q = 1).

Inverse regression step. The basic principle of the SIR method is to reverse the roles of y and x, that is, instead of regressing the univariate variable y on the multivariate variable x, the covariable x is regressed on the response variable y.

Let T denote a monotone (but not necessarily strictly monotone) transformation of y. Let µ = E(x) and Σ = V(x). We assumed that E((x ′ x) 2 ) < ∞ and that Σ is invertible. Under model (4.9) and (LC), [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] established the following geometric property: the centered inverse regression curve, E(x|T (y))µ as y varies, is contained in the linear subspace of R p spanned by ΣB. A straightforward consequence is that the covariance matrix

Γ := V(E(x|T (y)))
is degenerated in any direction Σ-orthogonal to Span(B). Therefore, the eigenvectors associated with the K non-null eigenvalues of Σ -1 Γ are some EDR directions.

Slicing step. [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] proposed a transformation T , called "slicing", which categorizes the response y into a new response with H > K levels. The support of y is partitioned into H non-overlapping slices s 1 , . . . , s h , . . . , s H . With such a transformation T , the subspace recovered through a slicing (based on the inverse x|T (y) function) may fall short of the space recovered through y in its entirety (based on the x|y function). However the main advantage of the slicing is that the matrix of interest is now written as

Γ = H h=1 p h (m h -µ)(m h -µ) ′
where p h = P(y ∈ s h ) and m h = E(x|y ∈ s h ).

Estimation process. In the usual statistical framework, when a sample {(x i , y i ), i = 1, . . . , n} is available, it is straightforward to estimate the matrices Σ and Γ, by substituting empirical versions of the moments for their theoretical counterparts, and therefore to obtain the estimation of the EDR directions. [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] showed that each of these estimated EDR directions converges to an EDR direction at rate √ n. Asymptotic normality of these estimated EDR directions has been obtained by [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF].

From a practical point of view, the choice of the slicing is discussed in [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF] or Saracco (2001). Since SIR theory makes no assumption about the slicing strategy, the user must choose the number H of slices and how to construct them. In practice, there are naturally two possibilities: to fix the width of the slices or to fix the number of observations per slice. This second option is often preferred, and from the sample point of view, the slices are often chosen such that the number of observations in each slice is as close to each other as possible. Note that H must be greater than K in order to avoid an artificial reduction of dimension and must be lower than ⌊n/2⌋ in order to have at least two observations in each slice (where ⌊a⌋ denotes the integer part of a). The choice of H is less sensitive than the choice of a smoothing parameter in nonparametric regression. This point is clearly illustrated in Liquet and Saracco (2012) with a graphical tool that allows the user to find simultaneously realistic values of the two parameters H and K, see Section 4. 3.3.5 for some details on this method.

The SIR estimates based on the first inverse moment have been studied extensively, see for instance [START_REF] Hsing | An asymptotic theory for sliced inverse regression[END_REF], [START_REF] Zhu | Asymptotics of sliced inverse regression[END_REF], Saracco (1999), [START_REF] Prendergast | Influence functions for sliced inverse regression[END_REF], [START_REF] Szretter | The sliced inverse regression algorithm as a maximum likelihood procedure[END_REF] among others for some asymptotic results. [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF] exhibited many features to popularize SIR. [START_REF] Chavent | A sliced inverse regression approach for a stratified population[END_REF] considered the case of a stratified population. In order to avoid the choice of a slicing in SIR, pooled slicing, kernel or spline versions of SIR have been investigated, see for example [START_REF] Zhu | Asymptotics for kernel estimate of sliced inverse regression[END_REF], [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF], [START_REF] Zhu | On spline approximation of sliced inverse regression[END_REF], [START_REF] Wu | Kernel sliced inverse regression with applications to classification[END_REF], Kuentz et al. (2010) or [START_REF] Azaïs | Optimal quantization applied to sliced inverse regression[END_REF]. However, these methods are hard to implement comparing to the basic SIR approach and are often computationally slow. Regularized versions for SIR have been proposed for high-dimensional covariates, see for instance [START_REF] Zhu | On sliced inverse regression with highdimensional covariates[END_REF], [START_REF] Scrucca | Class prediction and gene selection for DNA microarrays using regularized sliced inverse regression[END_REF], [START_REF] Li | Sliced inverse regression with regularizations[END_REF], Bernard-Michel et al. (2009b). Sparse SIR has been proposed by [START_REF] Li | Sparse sliced inverse regression[END_REF]. Hybrid methods of inverse regression-based algorithms have been also studied, see for example [START_REF] Gannoun | An asymptotic theory for SIR α method[END_REF] or Zhu et al. (2007).

Multivariate SIR

In the multivariate framework (that is when y ∈ R q with q > 1), [START_REF] Aragon | A gauss implementation of multivariate sliced inverse regression[END_REF], [START_REF] Li | Dimension reduction for multivariate response data[END_REF] considered several estimation methods of the EDR space based on SIR. Note that [START_REF] Barreda | Some extensions of multivariate SIR[END_REF] proposed extensions of the following multivariate SIR methods based on SIR α approach instead of SIR, where SIR α is a generalization of SIR which combines information from the first two conditional moments of x given T (y).

Complete slicing and marginal slicing approaches. In the complete slicing method, the SIR procedure is directly applied on y. To build slices of nearly equal sizes, the following recursive approach is used. The first component of y is sliced. Then, each slice is separately sliced again according to the next component of y, and so on. This extension of univariate SIR to multivariate y appears straightforward and the theoretical development can be formally carried over. Computation of such estimators suffers from the so-called curse of dimensionality when the dimension q of y is large (q ≥ 4). Note that [START_REF] Hsing | Nearest neighbor inverse regression[END_REF] proposed a version of SIR in which the slices are determined by the nearest neighbors approach and showed that the EDR directions can be estimated with rate √ n under general conditions. Moreover, [START_REF] Setodji | K-means inverse regression[END_REF] extended that univariate SIR to multivariate framework by introducing a new way of performing the slicing of y based on k-means method. The corresponding method is called k-means inverse regression (KIR).

A natural way to circumvent the curse of dimensionality of the complete slicing approach is proposed in the marginal slicing procedure which consists in applying the SIR method on a transformation of y depending on one's interest. For instance, it can be the mean or the median of the y (j) 's. One can also take the first few significant components of a principal component analysis of the y (j) 's to construct the slices. However, slicing a lower dimensional projection of y may not lead to recover as many EDR directions as slicing the entire y.

For these reasons, these two multivariate approaches (complete slicing and marginal slicing) are not completely satisfactory.

Pooled marginal slicing approach. The idea of the pooled marginal slicing (PMS) method is to consider the q univariate marginal SIR of each component y (j) of y on x and to combine the corresponding matrices of interest Γ (j) := V(E(x|T j (y (j) ))) in the following pooling: (4.11) for positive weights w j . It has been shown that the eigenvectors associated with the nonnull K eigenvalues of Σ -1 Γ P are EDR directions. [START_REF] Aragon | A gauss implementation of multivariate sliced inverse regression[END_REF] proposes to use two kinds of weighting for the w j 's: equal weights or weights proportional to the major eigenvalues found by a preliminary univariate SIR analysis of each component of y. Saracco (2005) obtained the asymptotic normality of the pooled marginal slicing estimator based on SIR α . [START_REF] Lue | Sliced inverse regression for multivariate response regression[END_REF] derived the asymptotic weighted chi-squared test for dimension. For j = 1, . . . , q, rather than constructing Γ (j) from y (j) , one can also build it from a linear combination τ ′ y of y. This method which is called projective resampling was introduced by Li et al. (2008). To ensure good performances, the number of linear combinations to handle should be greater than the sample size n.

Γ P = q j=1 w j Γ (j) ,
Some other multivariate SIR approaches. [START_REF] Bura | Estimating the structural dimension of regressions via parametric inverse regression[END_REF] introduced the parametric inverse regression that may easily adapt to multivariate response framework. [START_REF] Yin | Moment-based dimension reduction for multivariate response regression[END_REF] proposed a moment-based dimension reduction approach in this context. Moreover, in order to solve the dimensionality problem when p is large and to rationalize the slicing step, [START_REF] Li | Dimension reduction for multivariate response data[END_REF] presented an algorithm based on a duality between SIR variates and MP (most predictable) variates. The term "variate" denotes any linear combination of either the regressor x or the response variable y. The SIR variates are the variables b ′ x formed by an EDR direction b obtained with SIR. The MP variates θ ′ y are defined as those minimizing the ratio E[V(θ ′ y|x)]/V(θ ′ y), where V(θ ′ y|x) is the associated prediction mean squared error of the best nonlinear prediction E[θ ′ y|x] for the squared error loss. Equivalently, due to ANOVA identity, the MP variates can be found by maximizing the ratio V(E[θ ′ y|x])/V(θ ′ y), which conducts to the same eigenvalue decomposition as the SIR approach except for the exchanged roles of x and y. This twin relationship between SIR variates and MP variates underlies the development of the alternating SIR algorithm. The idea of the algorithm is to alternate computations of either θ or b respectively obtained by the slicing of SIR variates or MP variates constructed at the previous step. [START_REF] Li | Dimension reduction for multivariate response data[END_REF] proposed an iterative procedure for the alternating SIR and showed that choosing the canonical directions as an initial projection of the y's guarantees the convergence of the corresponding algorithm in a finite number of steps (equal to K, the number of EDR directions).

Sample version

Let {(x i , y i ), i = 1, . . . , n} be a sample of independent observations from model (4.9). Each y i is a q-dimensional random variable. We assume that the sample size n is larger than the dimension p of each covariate x i . Let x := 1 n n i=1 x i and Σ := 1 n n i=1 (x ix)(x ix) ′ be the empirical mean and the covariance matrix of the x i 's.

In order to estimate the matrix B, we have to estimate each p × K matrix B (j) with usual univariate SIR from the subsample {(x i , y (j) i ), i = 1, . . . , n} where y (j) i stands for the jth component of y i . To this end, let us assume that the support of y (j) is partitioned into a fixed number of slices denoted by s

(j) 1 , . . . , s (j) h , . . . , s (j) H (j) . Let p (j) h := P y (j) ∈ s (j) h and m (j) h := E x|y (j) ∈ s (j) h
. Thus, the matrix Γ (j) := H (j) h=1 p

(j) h (m (j) h -µ)(m (j) h -µ) ′ is estimated by Γ (j) := H (j) h=1 p(j) h ( m (j) h -x)( m (j) h -x) ′ with p(j) h := 1 n n i=1 I[y (j) i ∈ s (j)
h ] and m (j)

h := 1 np (j) h n i=1 x i I[y (j) i ∈ s (j)
h ] where I[.] is the indicator function. Then, the matrix B (j) is estimated by

B (j) := b(j) 1 , . . . , b(j) K ,
where the vectors b(j) k , k = 1, . . . , K are the Σ-orthonormal eigenvectors associated with the K largest eigenvalues of the matrix Σ -1 Γ (j) . It follows that the matrix B is directly estimated by

B := B (1) , . . . , B (q) .
Finally, a Σ-orthonormal estimated basis of the EDR space is given by the vectors v1 , . . . , vK defined as the eigenvectors associated with the K largest eigenvalues of the matrix B B ′ Σ and we write V := [v 1 , . . . , vK ].

Remark 9. Similarly to the pooled marginal slicing (PMS) presented in Section 4.3.2.2, MSIR relies on the univariate version of SIR, applied to each component of y. While both methods need estimates of Γ (1) , . . . , Γ (q) , estimates of B (1) , . . . , B (q) are only required by MSIR. Computing such estimates is useful to explore the relations between components of y as explained in Section 4.3.4.

An asymptotic result

The following assumptions are necessary to state our asymptotic result. Let n (j)

h := np (j) h
be the number of observations in the slice s (j) h .

• (A1) Observations {(x i , y i ), i = 1, . . . , n} are independently drawn from a given regression model.

• (A2) For each component y (j) of y, the support is partitioned into a fixed number H (j) of slices such that p (j)

h > 0 for h = 1, . . . , H (j) .

• (A3) For j = 1, . . . , q and h = 1, . . . , H (j) , n

h → ∞ as n → ∞.

Theorem 10. Under model (4.9) and assumptions (LC) and ( A1)-(A3), we have, for k = 1, . . . , K,

vk = v k + O p (n -1/2 ),
which means that for all ǫ > 0 there exist ζ ǫ > 0 and ñǫ ∈ N * such that for all n ≥ ñǫ

P √ n |v k -v k | > ζ ǫ < ǫ,
and then the estimated EDR space Span( V) converges in probability to the EDR space.

The poof is given in Appendix A.3.3.

Remark 10. Using Delta-method and asymptotic results of [START_REF] Tyler | Asymptotic inference for eigenvectors[END_REF] and [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF], it is possible to obtain the asymptotic normality of

√ n vec( B B ′ Σ) -vec(BB ′ Σ) ,
where vec(M) stands for the "vec" operator applied to matrix M. More precisely, this operator rearranges the p 2 elements of M in the form of a p 2 -dimensional column vector by stacking the p columns of M one under the other. Then, the asymptotic normality of the eigenprojector onto the estimated EDR space can be derived, as well as the asymptotic distribution of the estimated EDR directions vk , associated with eigenvalues assumed to be different (that is

λ 1 > • • • > λ K > 0).

A weighted version of MSIR

Following the idea used in pooled marginal slicing approach in which the matrix of interest Γ P is a weighted average of the marginal matrices Γ (j) , we can consider a weighted version of the multivariate SIR method introduced in this paper, named wMSIR hereafter. As it has already been proposed by [START_REF] Aragon | A gauss implementation of multivariate sliced inverse regression[END_REF] or [START_REF] Lue | Sliced inverse regression for multivariate response regression[END_REF], we shall use weights based on the proportion of eigenvalues corresponding to significant eigenvectors (which are EDR directions) in each marginal SIR (i.e. univariate SIR on each marginal component of y).

More precisely, for j = 1, . . . , q, let λ

k , k = 1, . . . , p be the eigenvalues of the eigendecomposition

Σ -1 Γ (j) v (j) k = λ (j) k v (j) k where λ (j) 1 ≥ λ (j) 2 ≥ • • • ≥ λ (j)
p . Let us define, for each component y (j) of y, the proportion of eigenvalues corresponding to significant eigenvectors:

π (j) = K k=1 λ (j) k p k=1 λ (j) k
. Let us also define π ⋆ = q j=1 π (j) . Then, we can introduce the following qK × qK matrix of weights

W = diag(W (1) , . . . , W (j) , . . . , W (q) ) with W (j) = π (j)
π⋆ I K for j = 1, . . . , q. Note that, from a theoretical point of view, under model (4.9) and (LC), the matrix of weights is such that W = 1 q I qK since λ (j) k = 0 for j = 1, . . . , q and k = K + 1, . . . , p.

The population version of wMSIR consists in noticing that the eigenvectors ṽ1 , . . . , ṽK associated with the K largest eigenvalues of the Σ-symmetric matrix BWB ′ Σ span the EDR space. We write V := [ṽ 1 , . . . , ṽK ]. To show this result, one can proceed analogously to the proof of Theorem 9.

Let λ(j) k be the kth eigenvalue of Σ -1 Γ (j) . The sample version of wMSIR is obtained by substituting the empirical matrices B, W and Σ for their theoretical counterparts B, W and Σ, where W = diag( W (1) , . . . , W (q) ) with, for j = 1, . . . , q, W (j) = π(j)

π⋆ I K , π(j) = K k=1 λ(j) k p k=1 λ(j)
k and π⋆ = q j=1 π(j) . Therefore, one can get the corresponding estimated EDR directions V := ṽ1 , . . . , ṽK .

Theorem 11. Under model (4.9) and assumptions (LC) and ( A1)-(A3), we have, for k = 1, . . . , K,

ṽk = ṽk + O p (n -1/2 ),
that is, the estimated EDR space Span( V) converges in probability to the EDR space.

The proof is given in Appendix A.3.4.

Discussion on the choice of K and H

Up to now, the dimension K of the EDR space was assumed to be known. However, in most applications based on real datasets, the number K of indices β ′ k x in model (4.7) is a priori unknown and hence must be determined from the data. In addition, the number of slices H in MSIR has to be chosen.

Several approaches to determine K have been proposed in the literature for univariate SIR. Some of them are based on hypothesis tests on the nullity of the last (p -K) eigenvalues, see for instance [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], [START_REF] Schott | Determining the dimensionality in sliced inverse regression[END_REF], [START_REF] Bai | A chi-square test for dimensionality for non-gaussian data[END_REF], [START_REF] Barrios | A bootstrap method for assessing the dimension of a general regression problem[END_REF] or [START_REF] Nkiet | Consistent estimation of the dimensionality in sliced inverse regression[END_REF]. In the multivariate response framework, [START_REF] Lue | Sliced inverse regression for multivariate response regression[END_REF] derived an asymptotic weighted chi-squared test for dimension adapted to the pooled marginal slicing estimator. In our case, a crude choice of the dimension can be also made by a visual inspection of the eigenvalues scree plot of the matrix B B ′ Σ: the idea is to determine the number of the significantly non-null eigenvalues.

In the univariate response model, Liquet and Saracco (2012) proposed to consider a risk function which can be replaced in this multivariate context by:

R h,k := E r(V k , Σ, V k , Σ) , (4.14)
where . . . , vk ] and h is the number of slices used to obtain V k and V k . This risk function only makes sense for any dimension k lower than or equal to the true dimension K of the EDR space. For the true dimension K, R h,K converges to one as n tends to infinity. For a fixed n, a reasonable way to assess whether an EDR direction is available is to graphically evaluate how much R h,k departs from one. From a computational point of view, consistent estimates of R h,k are required. Liquet and Saracco (2012) use a bootstrap estimator R h,k of this criterion in order to determine the pair (H, K) of parameters. The proposed graphical method consists in evaluating the R h,k values for all k = 1, . . . , p and some reasonable values of h, and in observing how much the criterion departs from one. The best choice will be the pair ( H, K) which gives a value of R h,k close to one, such that K << p in order to get an effective dimension reduction.

V k := [v 1 , . . . , v k ], V k := [v 1 ,
In practice, there is no objective criterion to find a trade-off between a large value of the criterion R h,k and a small value of the dimension K. Then, a visual expertise of the 3D-plot of the R h,k versus (h, k) allows the selection of the best value. It is also useful to provide, for each (h, k), the boxplots of the bootstrap replication of the squared trace correlation to have a look on how much the corresponding k-dimensional linear subspace is stable. Although boxplots of R h,k are also useful to determine the optimal number of slices H, wMSIR is not really sensitive to this parameter as shown in Section 4.3.5.

Analyzing components of y through MSIR

Recall that the estimate V (resp. V) of MSIR (resp. wMSIR) is computed from the estimated EDR directions B(j) associated with each component of y. From these estimates, it is straightforward to calculate the proximity measure rj := r( B(j) , Σ, V, Σ) (resp. rj := r( B(j) , Σ, V, Σ)) between each estimated marginal EDR space and the estimated common one, for j = 1, . . . , q. Then it is easy to sort these measures in descending order and to draw the associated scree plot. For j = 1, . . . , q, assuming model (4.9) and observing a low value of rj or rj could indicate an unprecise estimate of B (j) since rj and rj tends to 1 in probability as n goes to ∞. One can then withdraw the component y (j) of y to improve the accuracy of V or V.

In addition, assuming model (4.9) with a low dimensional common EDR space for the whole components of y does not always seem realistic in real data analysis. Therefore, applying any multivariate SIR method on y should not provide a suitable dimension reduction. However, it makes sense to assume that only groups of components of y rely on model (4.9) with small values of K, as in model (4.10). For this model, we propose a methodology in order to identify the variables y (j) which share the same EDR space. Thus, we obtain clusters of components on which applying a multivariate SIR approach is sensible. Note that performing marginal univariate SIR on each component y (j) of y leads to consistent estimates of each K-dimensional EDR space, since for l = 1, . . . , L, it is assumed that the rank of B l is equal to K. Recalling notations of Section 4.3.3.2, we obtain q 1 estimates Span( B(1) ), . . . , Span( B(q 1 ) ) of Span(B 1 ), q 2 estimates Span( B(q 1 +1) ), . . . , Span( B(q 1 +q 2 ) ) of Span(B 2 ), and so on, but values of q 1 , . . . , q L are unknown, as well as the number L of clusters.

For (j, j ⋆ ) ∈ {1, . . . , q} 2 , we define rj,j ⋆ := r( B(j) , Σ, B(j ⋆ ) , Σ). Without any loss of generality, assume that L = 2. Let us define (j 1 , j 2 , j 3 , j 4 ) ∈ {1, . . . , q 1 } 2 × {q 1 + 1, . . . , q 1 + q 2 } 2 .

We thus have the following Lemma.

Lemma 12. Under model (4.10) and assumptions (LC) and (A1)-(A3), rj 1 ,j 2 and rj 3 ,j 4 tend to 1 in probability.

The proof of this Lemma is given in Appendix A.3.5. Let us remark that, however, rj 1 ,j 3 does not converge to 1 in probability since r(B 1 , Σ, B 2 , Σ) < 1. This leads to the following criterion to cluster components of y: for (j, j ⋆ ) ∈ {2, . . . , q} × {1, . . . , j -1}, components y (j) and y (j ⋆ ) are put in the same cluster if rj,j ⋆ is close to 1. To do so, we can perform, for instance, a hierarchical ascending classification on the q × q (symmetrtic) matrix of the proximity measures rj,j ⋆ (with rj,j = 1 for j = 1, . . . , q). Other clustering procedures can be applied on this matrix, such as the multidimensional scaling together with the k-means method. In Sections 4.3.5-4.3.6, we give illustrations of a clustering step for simulated and real datasets, which clearly improves the estimation of the corresponding EDR spaces.

Remark 11. Computing the common estimated EDR space for each obtained cluster is not time consuming since B (1) , . . . , B (q) have already been computed. This is not the case for the k-means inverse regression (KIR) which requires new computations. In addition, applying PMS instead of MSIR or wMSIR on a cluster of y requires to store the p × p matrices Γ

(1) , . . . , Γ q) and summing some of them, which represent more computational time and more memory space than required by MSIR or wMSIR method, especially when p is large.

A simulation study

This section illustrates the ability of the proposed MSIR and wMSIR approaches, together with the diagnostic procedures on the components of y, to properly estimate EDR spaces.

The two following subsections respectively correspond to models (4.9) and (4.10).

Single EDR space model

Two simulation models are considered here. For a given sample size n and a dimension p, 100 replications of the covariate x are generated from the p-dimensional normal distribution N p (µ, Σ) with the same pair (µ, Σ). Here, µ is randomly generated from the N p (0, I p ) distribution and Σ = LL ′ + 0.1I p where L is a p × p matrix made of entries independently generated from the standard normal distribution N 1 (0, 1). For every model, for all j ≥ 1, it is assumed that ε (j) ∼ N 1 (0, 1). The dimension K is also assumed to be known in this section. We first study a special case of model (4.9) with K = 1. Secondly, MSIR and wMSIR are evaluated with a multiple indices model (K = 2).

Single index model. Consider the following single index model:

   y (1) = x ′ β 1 + ε (1) , y (2) = (x ′ β 1 ) 3 + 3ε (2) , y (3) = x ′ β 1 (1 + ε (3) ), (4.15)
with

β 1 = [β 1,1 , . . . , β 1,p ] ′ .
We choose for all i = 1, . . . , p, β 1,i = i I(i ≤ 5) + I(i > 5). We generate samples of size n = 100 from (4.15), with p = 20. Then, the EDR direction is estimated using the following methods, with H = 10 slices:

• univariate SIR for each component of y which produces estimates B (1) , B (2) and B (3) ,

• MSIR which gives the estimate V,

• wMSIR that leads to the estimate V,

• k-means inverse regression which provides the estimate V KIR ,

• pooled marginal slicing, leading to the estimate V PMS .

For each estimator B ∈ B (1) , B (2) , B (3) , V, V, V KIR , V PMS , we compute the squared trace correlation r( B) := r( B, Σ, B, Σ) between the estimated EDR space and the true EDR space. The closer to one is r( B), the better is the estimate. Note that for K = 1, the criterion r( B) corresponds to the squared cosine of the angle between B and β 1 . Boxplots of this criterion are drawn on Figure 4.5(a). It appears that B (3) exhibits low squared trace correlation. This phenomenon can be explained by the heteroscedasticity in the third marginal model of (4.15). Even if B (3) is necessary to compute V and V, the poor estimates of B (3) do not imply a significant loss in the squared trace correlations related to V and to V. One can also observe that the weighting in wMSIR seems to improve the estimation of the EDR space since values of r( V) are globally greater than those of r( V). This trend is confirmed in Figure 4.5(b) where the values of r( V) are plotted versus those of r( V). It appears that wMSIR seems to be uniformly better than MSIR in this simulation. In addition, Figure 4.5(a) shows that the pooled marginal slicing produces slightly better estimates than wMSIR in this example and that wMSIR outperforms the k-means inverse regression. The poor quality of the estimate B (3) can be observed directly from the simulated data. In Figure 4.6(a), we draw boxplots of values of rj for j = 1, 2, 3. Considering that the first quartile of the third boxplot is equal to 0.80 and that the minimum value of the others is equal to 0.96, it makes sense to withdraw the third component of y from the analysis for at least a forth of the datasets. Let V ⋆ be the wMSIR estimate built only from y (1) and y (2) . In Figure 4.6(b), we compare quality measures for V ⋆ and for the PMS estimate. We observe that the selection based on the rj 's improves performances of the wMSIR method, so that it produces better quality measures than the PMS without this selection step.

To study the sensitivity of wMSIR with respect to n and p, several samples are generated for various values of these parameters. In Figure 4.7(a), boxplots of 100 values of r( V) are drawn for each pair (n, p). Not surprisingly, estimated EDR spaces become closer to the true one when n increases. Moreover, one can also observe that estimates of the EDR space are more precise when p is small. This phenomenon can be explained by the fact that p × p matrices have to be estimated. Notice that for every value of (n, p) in Figure 4.7(a), estimates are very precise since most of the values of r( V) are greater than 0.85.

In previous analyses, the number of slices H was set to 10 to have enough slices to study the functions that link x ′ β 1 to each component of y and enough points in each slice. In Figure 4.8(a), we observe that for wMSIR, one can arbitrarily choose a number of slices between 8 and 38 and obtain an estimate of the EDR space which is as reliable as the one computed with H = 10.

Multiple indices model. Consider now a more complex model than model (4.15). It is defined by:

y (1) = exp(x ′ β 1 ) × (x ′ β 2 ) + ε (1) , y (2) = (x ′ β 1 ) × exp(x ′ β 2 ) + ε (2) , (4.16)
where

β 1 = [β 1,1 , . . . , β 1,p ] ′ , β 2 = [β 2,1 , . . . , β 2,p ] ′ , with β 1,i = i I(i ≤ 5) + I(i > 5) and β 2,i = (-1) i-1 (1 + I(i ∈ {3, 4})).
Note that, in model (4.16), we clearly have K = 2. Samples are generated from model (4.16) for various n and p. Then, these samples are used to estimate the corresponding EDR space with wMSIR with H = 10 slices. This leads to boxplots of r( V) displayed in Figure 4.7(b). We observe lower and more scattered values of r( V) than in Figure 4.7(a). The complexity of the link function between y and x and the greater dimension K are believable reasons for this phenomenon. Apart from this feature, Figure 4.7(b) provides identical evolutions of r( V) with n and p to those observed for the single index model in Figure 4.7(a). Note that for both model (4.15) and model (4.16), the behavior of MSIR and wMSIR when n and p vary are similar. That is why only results concerning wMSIR are drawn in Figure 4.7.

Examining Figure 4.8(b), it seems that the quality of wMSIR estimates is less connected to the chosen number of slices for model (4.15) than for model (4.16). For the latter, performances of wMSIR are nevertheless similar for values of H from 6 to 12.

Multiple EDR spaces model

Consider model (4.10) with q = 12, p = 20 and K = 1. Define for i = 1, . . . , p, the ith component of respectively β 1 , β 2 and β 3 as β 1,i := i I(i ≤ 5) + I(i > 5), β 2,i := 6i + 5 i 5 and β 3,i := (-1) i-1 (1 + I(i ∈ {3, 4})). The random variable y is drawn from the following model: for the considered values of j, y (3j) is the third component of each model (4.15) embedded in model (4.17), these inequalities can be explained by the heteroscedasticity in this component which leads to imprecise estimates of the relevant EDR direction as pointed out on Figure 4.5(a). In addition, on Figure 4.9, component y (11) could be unexpectedly grouped with y (4) and y (5) or with y (8) .

In this example, the interpretation of Figure 4.9 can easily be done because components of y are already clustered in the definition of the model. In other words, every component between two others that are related to the same EDR space belongs to a marginal model based on this EDR space. In practical cases, the components of y may not be ordered that way which means that the corresponding representation of the squared trace correlations may be cluttered. To tackle this problem, we use an agglomerative hierarchical clustering algorithm based on the dissimilarity 1rj,j ⋆ between B(j) and B(j ⋆ ) .

A model with 5 clusters. We generate a new sample of size n = 1000 from model (4.17) with θ 1 = 1 and θ 2 = 0 which leads to l = 5 clusters: y (1) , y (2) , y (3) , y (7) , y (8) , y (9) , y (4) , y (5) , y (6) , y (10) , y (11) , y (12) .

The hierarchical clustering algorithm applied to the estimates B (j) for j = 1, . . . , q and produces the dendrogram of Figure 4.10(a). A classification directly based on this procedure would not be really accurate. For instance, in model (4.17), y (11) and y (10) belong to two different clusters. On Figure 4.10(a), to divide y (11) and y (10) into two groups, the tree has to be cut at a level around 0.75. This implies grouping y (1) , y (7) , y (2) and y (8) together and putting y (3) and y (9) in another group while components of both groups are actually related with the same EDR space. However, the dendrogram of Figure 4.10(a) allows to order components of y in such a way that those which are linked by high squared trace correlation are close from each other. Thus, in Figure 4.10(b), we displayed values of rj,j ⋆ for j ⋆ ≺ j where ≺ denotes the ordering in Figure 4.10(a). It becomes clear, then, that components y (1) , y (7) , y (2) , y (8) , y (3) and y (9) should be clustered in the same group. Note that another cluster containing y (6) , y (4) and y (5) can be made from Figure 4.10(b). Not surprisingly, we observe again that the squared trace correlation between components y (3) , y (6) , y (9) and the components that correspond to the group they respectively belong to is quite low. Finally, components y (12) , y (11) and y (10) appears to form three distinct clusters.

Real data illustrations 4.3.6.1 Remote sensing data

As an illustration, we consider a nonlinear inverse problem in remote sensing. The goal is to estimate the physical properties of surface materials on the planet Mars from hyperspectral data. The method is based on the estimation of the functional relationship between some physical parameters x and observed spectra y. We refer to Bernard-Michel et al. (2009a) for further details. We focus on an observation of the south pole of Mars at the end of summer 2003, collected by the French imaging spectrometer OMEGA on board the Mars Express Mission. A detailed analysis of this image [START_REF] Douté | South pole of Mars: Nature and composition of the icy terrains from Mars Express OMEGA observations[END_REF]) revealed that this portion of Mars mainly contains water ice, carbon dioxide and dust. This led to the physical modeling of individual spectra with a surface reflectance model y = g(x). The p = 3 parameters x (1) , x (2) and x (3) are respectively the proportion of carbon dioxide, the proportion of dust, and the grain size of water ice. Let us note that the proportion of water is equal to 1x (1)x (2) . Each spectra y is made of q = 352 wavelengths. The link function g has no close-form expression, but it can be computed thanks to a dedicated software. This yields the simulation of a sample {(x i , y i ), i = 1, . . . , n} of size n = 6400.

We ran the clustering procedure described in the above section associated to wMSIR with H = 10 slices. The clustering results are depicted on Figure 4.11. Two clusters of wavelengths have been identified, corresponding to two different orientations of ṽ1 . It appears on Figure 4.12(a) that, in the first cluster, only the proportion of dust x (2) has an important contribution to the EDR direction. In the second cluster, the estimated EDR direction is close to (1, 1, 0). This corresponds to the index x (1) + x (2) which is equal to one minus the proportion of water. Let us note that the grain size of water ice x (3) does not appear in the EDR directions. Figure 4.12(b) permits to visualize the clustering of the wavelengths: It reveals which wavelengths are more sensible to the presence of water ice or dust.

Minneapolis elementary schools data

Another dataset which is widely studied in the dimension reduction context with a multivariate response is related to test results of students in Minneapolis elementary schools. These data are for example presented in [START_REF] Shao | Partial central subspace and sliced average variance estimation[END_REF]; [START_REF] Cook | A model-free test for reduced rank in multivariate regression[END_REF]; [START_REF] Yin | Moment-based dimension reduction for multivariate response regression[END_REF]. The response variable y is made of q = 4 components. The first (resp. third) component is the proportion of pupils scoring below the average on a fourth (resp. sixth) grade test. The second (resp. fourth) one is the proportion of marks above the average. Following [START_REF] Yin | Moment-based dimension reduction for multivariate response regression[END_REF], we would like to explain y with a 8-dimensional variable x. The seven first components x are called x (1) , . . . , x (7) and are respectively the squared root of the percentages of children receiving an aid called AFDC, children who do not live with both parents, people in the area of a school who completed high school, people who suffer for poverty, minority, mobility and pupils who attend school regularly. The last component of x, named x (8) , is the mean number of pupils for each teacher. Note that [START_REF] Yoo | Iterative optimal sufficient dimension reduction for conditional mean in multivariate regression[END_REF] analyzed the variables x (1) , x (2) and x (3) while [START_REF] Shao | Partial central subspace and sliced average variance estimation[END_REF] focused on some increasing functions of these variables and [START_REF] Cook | A model-free test for reduced rank in multivariate regression[END_REF] considered

x (1) , x (2) , x (3) , x (4) and x (8) . We first compute V with p = 8 components of x. We then use a stepwise algorithm with the AIC to perform a linear regression of x ′ V on x and sort the components of x with respect to how much they explain x ′ V. The first three sorted variables are x (1) , x (2) , and x (3) . Besides, regressing x ′ V on x (1) , x (2) , x (3) produces an adjusted coefficient of determination of 95.42%. This encourages us to take x = (x (1) , x (2) , x (3) ) ′ .

Working with this 3-dimensional covariate and with H = 8, we compute B (j) for j = 1, . . . , 4, with K = 1, which is the size of the dimension chosen in [START_REF] Shao | Partial central subspace and sliced average variance estimation[END_REF], [START_REF] Cook | A model-free test for reduced rank in multivariate regression[END_REF] and [START_REF] Yin | Moment-based dimension reduction for multivariate response regression[END_REF]. We then construct the matrix B B ′ Σ. Its eigenvalues are 0.96, 0.04 and 0.01 which confirms the choice of K = 1. We can not build several groups of components of y since the least value of r(j, j ⋆ ) for (j, j ⋆ ) ∈ {1, . . . , 4} is equal to 0.75, neither we can withdraw a component from y since the least value of rj for j ∈ {1, . . . , 4} is equal to 0.92. In addition, we find V = (0.673, -0.406, -0.528) ′ . Let V Y be the EDR direction found by [START_REF] Yoo | Iterative optimal sufficient dimension reduction for conditional mean in multivariate regression[END_REF]. We have that r( V, Σ, V Y , Σ) = 0.97.

The signs of the elements of V make also sense compared to results from [START_REF] Shao | Partial central subspace and sliced average variance estimation[END_REF] and [START_REF] Cook | A model-free test for reduced rank in multivariate regression[END_REF].

Concluding remarks

In this paper, we proposed the new multivariate SIR approaches MSIR and wMSIR for estimating the EDR space. The idea consists in performing first several marginal SIR analyzes. A common EDR space is then deduced from the marginal ones by maximizing the proximity criterion defined in (4.12). This optimization problem benefits from a closed-form solution.

Let us highlight that this two-step approach can be applied to any variant of the SIR method. For instance, it permits to build multivariate regularized SIR approaches from Bernard-Michel et al. (2009b), multivariate SIR α approaches from [START_REF] Gannoun | An asymptotic theory for SIR α method[END_REF] or multivariate kernel SIR methods from [START_REF] Wu | Kernel sliced inverse regression with applications to classification[END_REF]. To avoid the choice of the number H of slices, one may also consider building a multivariate version of the CUME procedure from Zhu et al. (2010a).

MSIR and wMSIR can also be run with a graphical procedure that cluster components of the response variable depending on the EDR space they are related with. Therefore, the approach we described allows to deal with datasets that come from a model that includes several different EDR spaces. It is then possible to estimate each of them from clusters of components of the response variable rather than blindly apply a multivariate SIR procedure on the whole variable. In addition, we noticed that estimating the EDR space does not cost a significant amount of computational time when it is done in the context of this clustering procedure. R codes are available on request from the authors.

Comparison of sliced inverse regression approaches for underdetermined cases

This section is an article which is accepted for publication in Journal de la Société Française de Statistique. It was written with Benoît Liquet and Jérôme Saracco. Additional contents were nonetheless added in the simulation study of Section 4.4.4.

Introduction

For a univariate response variable y and a multivariate covariate x ∈ R p , the semiparametric regression model

y = f (x ′ β 1 , . . . , x ′ β K , ε) (4.18)
is an attractive dimension-reduction approach to model the effect of the p-dimensional covariates x on y. Let µ = E(x) and Σ = V(x). The error term ε is assumed to be independent of x. Since the link function f (•) is an unknown smooth function, the parameters β k ∈ R p are not entirely identifiable, only the linear subspace spanned by the β k 's can be identified without additional assumptions. [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] called this subspace the effective dimension reduction (EDR) subspace. Moreover any direction belonging to this subspace is called an EDR direction. If the β k 's are assumed linearly independent, the EDR subspace is then a K-dimensional linear subspace of R p .

Other authors refer to this subspace as the dimension reduction subspace (DRS) or the central subspace (which is defined as the smallest DRS), see Cook (1998) for more details. When the dimension p of x is high and when we have little knowledge about the structure of the relationship between the response and the covariates, this semiparametric regression model is a nice alternative to parametric modeling (since it is really difficult to have knowledge about the structure of the relationship between the response and the covariates) and non-parametric modeling (which suffers from the well-known curse of dimensionality due to the data sparseness in the domain of x). The idea of dimension reduction in model (4.18) is intuitive because it aims at constructing a low dimensional projection of the covariate without losing information to predict the response y. If the dimension K of the EDR subspace is sufficiently small, it facilitates data visualization and explanation and it alleviates the curse of the dimensionality to non-parametrically estimate f with usual approaches such as kernel or splines smoothing (when the error term is additive).

In this semiparametric regression model (4.18), an important purpose is to estimate the EDR subspace from a sample {(x i , y i ), i = 1, . . . , n}. Most of the existing approaches are usually based on the eigendecomposition of a specific matrix of interest. The most popular one is the sliced inverse regression (SIR) introduced by [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] and [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] , respectively for single index models (K = 1) and multiple indices models (K ≥ 1). Among alternative methods there are SIR-II, see [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF]; [START_REF] Yin | Asymptotic distributions for dimension reduction in the sir-ii method[END_REF] for instance, and sliced average variance estimation (SAVE), see Zhu and Zhu (2007); [START_REF] Li | Asymptotics for sliced average variance estimation[END_REF] for example. These approaches require the inverse of Σ. Then, from a practical point of view, it is necessary to inverse an estimate Σ of Σ. Define xi = (x iμ) ∈ R p for i = 1, . . . , n, with μ = 1 n n i=1 x i . A usual (biased) variable x, the covariable x is regressed on the response variable y. The price we have to pay to succeed in inverting the role of x and y is an additional assumption on the distribution of x, named the linearity condition (described hereafter).

Usual SIR estimate is based on the first moment E(x|y). It has been initially introduced by [START_REF] Duan | Slicing regression: a link-free regression method[END_REF] for single index model and by [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] for multiple indices model. SIR approaches have been extensively studied, see for instance [START_REF] Carroll | Measurement error regression with unknown link: dimension reduction and data visualization[END_REF]; [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF]; Zhu et al. (2007); [START_REF] Bercu | A new approach of recursive and non recursive sir methods[END_REF]; [START_REF] Azaïs | Optimal quantization applied to sliced inverse regression[END_REF] among others.

Let us now recall the geometric property on which SIR is based. Let us introduce the linearity condition:

(LC) : ∀b ∈ R p , E(x ′ b|x ′ β 1 , . . . , x ′ β K ) is linear in x ′ β 1 , . . . , x ′ β K . (4.20)
Note that this condition is satisfied when x is elliptically distributed (for instance normally distributed). The reader can find an interesting discussion on this linearity condition in [START_REF] Chen | Can SIR be as popular as multiple linear regression?[END_REF].

Assuming model (4.18) and (LC), [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] showed that the centered inverse regression curve is contained in the linear subspace spanned by the K vectors Σβ 1 , . . . , Σβ K . Let T denote a monotonic transformation of y. He considered the eigendecomposition of the Σ-symmetric matrix Σ -1 M where M = V(E(x|T (y))). Straightforwardly the eigenvectors associated with the largest K eigenvalues of Σ -1 M are some EDR directions.

Slicing step. To easily estimate the matrix M , [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF] proposed a transformation T , called a slicing, which categorizes the response y into a new response with H > K levels (in order to avoid an artificial reduction of dimension). The support of y is partitioned into H non-overlapping slices s 1 , . . . , s h , . . . , s H . With such transformation T , the matrix of interest M can be now written as M = H h=1 p h (m h -µ)(m h -µ) ′ where p h = P(y ∈ s h ) and m h = E(x|y ∈ s h ).

Estimation process. When a sample {(x i , y i ), i = 1, . . . , n} is available, matrices Σ and M are estimated by substituting empirical versions of the moments for their theoretical counterparts. Let

M = H h=1 ph ( mh -μ)( mh -μ) ′ , (4.21) where ph = 1 n n i=1 I[y i ∈ s h ] and mh = 1 np h n i=1 x i I[y i ∈ s h ].
Therefore the estimated EDR directions are the eigenvectors associated with the K largest eigenvalues of Σ -1 M . They span the K-dimensional estimated EDR subspace The convergence at rate √ n and the asymptotic normality of estimated EDR directions have been obtained, see [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF]; [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF] for instance.

The choice of the slicing T is discussed in [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF]; [START_REF] Kötter | Smoothing and regression. Approaches, computation and application[END_REF]; Saracco (2001) but, theoretically, there is no optimal one. In practice, we fix the number of observations per slice to ⌊n/H⌋ where ⌊a⌋ stands for the integer part of a. If the sample size n is not proportional to the number H of slices, some slices will then contain ⌊n/H⌋ + 1 observations. Note that, in order to avoid the choice of a slicing, alternative SIR methods have been investigated. For instance, one can mention kernel-based methods of SIR proposed by [START_REF] Zhu | Asymptotics for kernel estimate of sliced inverse regression[END_REF] or [START_REF] Aragon | Sliced Inverse Regression (SIR): an appraisal of small sample alternatives to slicing[END_REF]. However, these methods are hard to implement and are computationally slow. Moreover, [START_REF] Bura | L 1 -Statistical Procedures and Related Topics, chapter Dimension reduction via parametric inverse regression[END_REF] and [START_REF] Bura | Estimating the structural dimension of regressions via parametric inverse regression[END_REF] proposed a parametric version of SIR.

Concerning the determination of the dimension K of the EDR subspace (which is unknown in practice), several works are available in the literature, see for example [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF]; [START_REF] Schott | Determining the dimensionality in sliced inverse regression[END_REF]; [START_REF] Ferré | Determining the dimension in sliced inverse regression and related methods[END_REF]; [START_REF] Bai | A chi-square test for dimensionality for non-gaussian data[END_REF]; [START_REF] Liquet | Application of the bootstrap approach to the choice of dimension and the α parameter in the SIR α method[END_REF] among others.

Standardized version for SIR. Another way to obtain a basis of the EDR subspace is to consider the eigendecomposition of Σ -1/2 M Σ -1/2 , that is the eigendecomposition of M * = V(E(z|T (y)) where z = Σ -1/2 (xµ) is the standardized version of the covariate x.

For the multiple indices model (4.18), we then focus on the first K eigenvectors η 1 , . . . , η K associated with the largest K eigenvalues of the I p -symmetric matrix M * . Transforming back to the original scale, the vectors Σ -1/2 η k , k = 1, . . . , K are in the EDR subspace. Their estimation procedure is a straightforward replication of the previous estimation process using M * = Σ -1/2 M Σ -1/2 .

Two existing SIR methods when n < p

As previously mentioned, the rank of Σ implies that this matrix is singular when n < p and ill-conditioned when n ≈ p. In this section we present two methods, respectively from [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] and [START_REF] Li | Sliced inverse regression with regularizations[END_REF], to tackle these cases.

RSIR: A modified estimated variance matrix. [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] introduce an upgrade of the SIR method, called RSIR, that relies on a modification of Σ such that the result can be inverted. This leads to the following estimate of Σ:

Σ(s) = Σ + sI p ,
where s is a positive real parameter and I p is the p × p identity matrix. For a given matrix A, let A 2 = Trace(A ′ A). To find a suitable s, [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] propose to minimize the mean squared error

L(s) = K k=1 Trace(V( βk (s))) + K k=1 E( βk (s)) -β k 2 ,
where β(s) = ( β1 (s), . . . , βK (s)) is the matrix of the K first generalized eigenvector of M and Σ(s), which is built such that the constraint βk (s) ′ Σ(s) βk (s) = I[k = k] is verified for all (k, k) ∈ {1, . . . , K} 2 . More details about generalized eigenvectors can be found in Section 4.4.2.3. Because β k is unknown, the authors replaced it with E( βk (s 0 )) in the expression of L(s), to obtain an approximation L(s). Note that the parameter s 0 has to be sufficiently small in order for E( βk (s 0 )) to be close to β k . In practice, s 0 is chosen equal to 0. Variances and expectations in L(s) are then estimated with bootstrap samples, which leads to an estimate L(s) of L(s). Remark that estimating E( βk (s 0 )) for s 0 = 0, implies using SIR with Σ. To do so, [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] apply the QZ algorithm (see Section 4.4.2.3 for details). The optimal regularization parameter is then given by s opt = arg min s L(s).

The corresponding matrix of estimated EDR directions is finally defined by βRSIR = β(s opt ).

SR-SIR: A ridge sliced inverse regression. We describe here the SR-SIR method from [START_REF] Li | Sliced inverse regression with regularizations[END_REF]. When Σ is invertible, let β = ( β1 , . . . , βK ) be the p × K matrix made of the eigenvectors of Σ -1 M . According to [START_REF] Li | Sliced inverse regression with regularizations[END_REF] (see also [START_REF] Setodji | K-means inverse regression[END_REF]), β also satisfies

β, V = arg min u,v H h=1 ph ( mh -μ) -Σuv h 2 , (4.22)
with v = (v 1 , . . . , v h ) and where the minimum is taken over the respective sets of p × K matrices and K × H matrices. Note that this equation is also defined when Σ is not invertible. From (4.22), the authors proposed thus a ridge version of the estimator β for a given regularization parameter s:

β(s), V (s) = arg min u,v G s (u, v). (4.23)
where

G s (u, v) = H h=1 ph ( mh -μ) -Σuv h 2 + s u 2
In practice, β(s) can be then obtained from (4.23) with an alternating least-squares algorithm even when n < p. The SR-SIR method rely on a generalized crossvalidation criterion to find the optimal regularization parameter s opt (see [START_REF] Li | Sliced inverse regression with regularizations[END_REF] for details). Finally the matrix of estimated EDR directions is defined by βSR-SIR = β(s opt ).

Remark 12. The exisence of a solution for (4.23) is not proved as explained by [START_REF] Bernard-Michel | A note on sliced inverse regression with regularizations[END_REF]. Indeed, assume that β(s), V (s) is such a solution and that β(s) is not the null vector, we then have

G s 1 2 β(s), 2 V (s) < G s β(s), V (s) ,
which contradicts the fact that β(s), V (s) verifies (4.23). This encourages [START_REF] Bernard-Michel | A note on sliced inverse regression with regularizations[END_REF] to replace (4.23) with the following optimization problem:

β(s), V (s) = arg min u,v H h=1 ph ( mh -μ) -Σuv h 2 + s uv W 1/2 2 ,
where W = diag(p 1 , . . . , pH ). The value of β(s) in this problem is actually the estimate of the RSIR method, for a regularization parameter s.

SIR-QZ: Solving the generalized eigenvalues problem in SIR

When Σ is regular, usual SIR estimates of the EDR directions are eigenvectors of Σ -1 M . This eigendecomposition is actually a special case of a generalized eigenvalues problem which consists in finding real numbers λ and non-null vectors v such that:

M v = λ Σv. (4.24)
The generalized Schur decomposition. When Σ is singular the generalized eigenvalue problem can still be solved if the function λ → Mλ Σ behave properly. We call this function a matrix pencil. In this section, we present the QZ algorithm which allows us to find couples (λ, v) that verify (4.24) for a wide range of matrix pencils including some with singular matrices Σ. The QZ algorithm can be viewed as an extension of the QR algorithm and was proposed by [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF]. The reader can refer to chapter 7 of [START_REF] Golub | Matrix computations[END_REF] for details. A brief description of this algorithm is provided in the following.

Notice that if we have two invertible matrices Q and Z, then finding λ and v in (4.24) is equivalent to find λ and w in

Q M Zw = λQ ΣZw, (4.25)
and to set v = Zw. Similarly to the QR algorithm that is designed to find the Schur decomposition of a matrix in order to compute its eigenvalues, the QZ algorithm aims at finding unitary matrices Q and Z such that Q M Z and Q ΣZ are upper triangular, for square matrices M and Σ. Such a transformation is called a generalized Schur decomposition. When working with complex matrices, Q and Z always exist (see [START_REF] Golub | Matrix computations[END_REF], Theorem 7.7-1). Possible values of λ that verify (4.24) are such that det( Mλ Σ) = 0, and such that det(Q( Mλ Σ)Z) = 0. The latter determinant is the product of the diagonal elements of Q( Mλ Σ)Z since it is an upper triangular matrix. Hence, the generalized eigenvalues of (4.24) are the ratios of the diagonal elements of Q M Z to the ones of Q ΣZ, provided that the diagonal elements of Q ΣZ are not equal to zero. More specifically, this can be seen with the following formula [START_REF] Golub | Matrix computations[END_REF], Theorem 7.7-1)

det( M -λ Σ) = det(Q ′ Z ′ ) p j=1 (t j -λu j ), (4.26)
where t 1 , . . . , t p and u 1 , . . . , u p are the respective diagonal elements of Q M Z and Q ΣZ. Notice that the generalized Schur decomposition only produces complex upper triangular matrices Q M Z and Q ΣZ. However, there is a similar available decomposition for real matrices M and Σ (see Appendix A.2.1 or [START_REF] Golub | Matrix computations[END_REF] for details).

Estimating the indices X ′ β using the QZ algorithm. Equation (4.26) implies that if it exists j ∈ {1, . . . , p} such that t j = u j = 0, then det( Mλ Σ) = 0 for all λ ∈ C, and trying to choose eigenvectors corresponding to the greatest eigenvalues to estimate the EDR directions does not make sense. Numerically, for any j ∈ {1, . . . , p}, due to rounding errors, t j and u j are almost always different from 0, but if both their absolute value are too small, det( Mλ Σ) is sufficiently unstable to call its value in question. As a consequence, every computed λ j can be wrong. For similar reasons, if |u j | is too small for a given j ∈ {1, . . . , p}, t j /u j should not be considered as an eigenvalue. These remarks and the regularization procedure in [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] lead to the algorithm we provide in Appendix A.2.2 to find an estimate βQZ of the EDR directions. Let X = (x 1 , . . . , x n ). When n is sufficiently smaller than p and H > K, a generalized eigenvector v of (4.24) is such that the n indices X ′ v only takes H distinct values, as explained in Appendix A.2.3. In practice, the regularization parameter s of the algorithm of Appendix A.2.2 is small and we can distinguish easily H clusters in the values of X ′ βQZ in Figure A.3 in Appendix A.2.3. This is a drawback of our approach when n < p since the values of X ′ β are a priori distinct. To circumvent this shortcoming, we compute several βQZ with different number of slices H 1 , . . . , H N H . Let denote the corresponding estimates βQZ,1 , . . . , βQZ,N H . We would like to find a K-dimensional subspace of R n which is as close to the KN H points of the matrix X ′ ( βQZ,1 , . . . , βQZ,N H ) as possible. Thus, we would choose a basis γ of this subspace as an estimate of X ′ β. This leads us to consider the following equation:

γ, δ = arg min γ,δ X ′ ( βQZ,1 , . . . , βQZ,N H ) -γδ 2 (4.27)
where the minimum is taken over the respective sets of n × K matrices and K × KN H matrices and each column of γδ is the approximation of the corresponding column of X ′ ( βQZ,1 , . . . , βQZ,N H ) in the K-dimensional subspace spanned by γ. A solution of (4.27) is given by a principal component analysis (see Besse (2012), p80-81). Note that their exists an infinite number of bases γ which span the optimal K-dimensional subspace. Thus, the solution provided by the principal component analysis is just one of them. We call the whole approach SIR-QZ.

Remark 13. When n is smaller enough than p, finding a satisfying estimate of β may not be possible. For example, if K = 1, n < p and if the columns of X are not linearly dependent, there are infinitely many solutions u of the system X ′ u = X ′ β for a given estimate β of β. Recalling the underlying model (4.18), there is no reason why β should be a better estimate of β than any of these solutions. That is why, when n < p, we focus on estimates of X ′ β rather than on β itself, in (4.27).

SIR-MP: a generalization of the inverse for singular matrices

We describe in this section a method which mimics the SIR approach developed for a functional covariate.

Dimension reduction in functional regression. In the functional SIR context, x is an explanatory functional variable (assumed square integrable in order to have its covariance operator well-defined) while y is still a real response variable. In this context, while the covariance operator of x is invertible, it has unbounded inverse so that its estimator is ill-conditioned. Then several methods have been proposed when the covariance operator does not need to be inverted.

One of them consists in using the eigendecomposition of M + Σ instead of Σ -1 M , where M + is the Moore-Penrose generalized inverse of M , also called Moore-Penrose pseudoinverse of M . In the particular context of functional sliced inverse regression, the reader can find a discussion on the fact that the eigenvectors of Σ -1 M are eigenvectors of M + Σ, in [START_REF] Ferré | Reply to the paper : "A note on smoothed functional inverse regression[END_REF] and references cited therein.

Let us now focus on an alternative approach introduced by Amato et al. (2006). They used the fact that Σ -1/2 M Σ -1/2 is a finite rank operator, where Σ (resp. M ) stands here for the covariance operator of x (resp. E(x|T (y))) in this functional context. The eigenvectors of this operator are eigenvectors of Σ 1/2 M + Σ 1/2 . The authors claimed that the reason of their approach is that a smooth estimate of M produces more stable estimates of the eigenvalue decomposition of M than that of the empirical estimate of Σ. Thus the eigenfunctions η 1 , . . . , η K associated with the smallest K eigenvalues α 1 , . . . , α K of Σ 1/2 M + Σ 1/2 are also the eigenfunctions associated with the largest eigenvalues of Σ -1/2 M Σ -1/2 equal to 1/α k for k = 1, . . . , K. In order to transform back to the original scale, we can not use the transformation Σ -1/2 η k . A basis of the (functional) EDR space is instead given by b k = M + Σ 1/2 η k for k = 1, . . . , K.

(4.28)

We provide in Appendix A.3.6 a brief proof of this result of [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF].

Adaptation for multivariate real covariates. In the context of our paper (that is, n < p < ∞), we will evaluate how the functional SIR procedure behaves in the multivariate framework. To do this, we simply substitute the operators of covariance M and Σ by the estimates M and Σ previously defined in (4.19) and (4.21). The resulting estimated directions are:

b k = M + Σ 1/2 η k for k = 1, . . . , K,
where the η k 's are the eigenvectors of Σ 1/2 M + Σ 1/2 associated with the smallest eigenvalues (among those not structurally equal to zero, see Remark 14 for details). This adaptation is called SIR-MP.

Remark 14. For a p-dimensional covariate, the p × p matrix Σ 1/2 M + Σ 1/2 is symmetric positive semidefinite and its rank r is at most equal to H -1 when H < n < p. Therefore, the eigenvalues of Σ 1/2 M + Σ 1/2 are such that α1 ≥ • • • ≥ αr > 0, and the geometric multiplicity of the eigenvalue zero is equal to p-r by construction. Thus we are interested in the eigenvectors η k associated with the K eigenvalues αr , . . . , αr-K+1 .

Selecting relevant components of x which are linked with y

Let β j,k denote the jth element of the EDR direction β k , for k = 1, . . . , K and j = 1, . . . , p. If β j,. = (β j,1 , . . . , β j,K ) is the null vector then the jth component of x does not have any effect on y. Finding such components is an important concern when n < p because it allows x ∈ R p to be reduced to x ⋆ ∈ R p ⋆ , where p ⋆ < p, without any loss of information. If in addition p ⋆ is less enough than n, the EDR directions can then be accurately estimated with a classical SIR procedure applied on y and x ⋆ . In Section 4.4.3.1, we describe the methods from [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] and [START_REF] Li | Sliced inverse regression with regularizations[END_REF] to determine which β j,. are null. In Section 4.4.3.2, we introduce another method to solve this problem based on proximity measures between models with only a few components of x and the initial model (in which every component of x is taken into account).

4.4.3.1 Review of existing selection procedures RSIR: Bootstrap estimates and a chi-squared test. Let βRSIR,j,k be the elements of the matrix βRSIR . [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] claim that for j = 1, . . . , p, the vector βRSIR,j,. = ( βRSIR,j,1 , . . . , βRSIR,j,K ) follows asymptotically a multivariate normal distribution with mean β j,. (s opt ) and covariance matrix Γ j . Provided that Γ j can be inverted, if β j,. (s opt ) is the null vector then β′ RSIR,j,. Γ -1 j βRSIR,j,. follows asymptotically a chi-squared distribution with K degrees of freedom. This encouraged [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] to use a chi-squared test on β′ RSIR,j,. Γ -1 j βRSIR,j,. to select which components of x have effect on y, where Γ j is an estimate of Γ j computed from bootstrap estimates of β j,. (s opt ).

Difficulties using this procedure could arise when n < p because the distribution under the null hypothesis of this test is asymptotic. In addition, it leads to inferences about the vectors β j,. (s opt ) for j = 1, . . . , p which are not necessary the same than for β j,. . SR-SIR: A Lasso method. From β(s opt ) and (ṽ 1 , . . . , ṽH ) = V (s opt ), given in equation (4.23), [START_REF] Li | Sliced inverse regression with regularizations[END_REF] propose to minimize the following expression under a constraint on the L 1 -norm of the vector ϕ:

G(ϕ) = H h=1 ph ( mh -μ) -Σdiag(ϕ) β(s opt )ṽ h 2 .
This leads to the following optimization problem for a parameter τ > 0:

φτ = arg min ϕ {G(ϕ)} , s.t. |ϕ| ≤ τ
where the minimum is taken over the set of vectors ϕ of length p. The Lasso procedure [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF]) can be used to find φτ . [START_REF] Li | Sliced inverse regression with regularizations[END_REF] consider that a component of x that corresponds to a zero in φτ does not have any effect on y. Let p τ be the number of non-null elements of φτ . In practice, choosing τ implies choosing the amount of selection provided by φτ . To do so, [START_REF] Li | Sliced inverse regression with regularizations[END_REF] propose to use classical model selection criteria. More specifically, this involves minimizing one of the following expression over a set of tested values of τ :

AIC = pH log G( φτ ) pH + 2p τ , BIC = pH log G( φτ ) pH + log(pH)p τ , RIC = (pH -p τ ) log G( φτ ) pH -p τ + p τ (log(pH) -1) + 4 pH -p τ -2 .

CSS based on SIR: Closest submodel selection for SIR methods

The idea of the procedure described here is to select submodels of (4.18) with only a given number p 0 of components of x which are the closest to the initial one. The components of x that appear the most in these submodels are asserted to have an effect on y.

Let Y = (y 1 , . . . , y n ) ′ . To do this, we propose the following algorithm.

Initialize

p 0 ∈]1, p[, N 0 ∈ N * and ζ ∈]0, 1[ or ρ ∈]0, 1[.
Step 1. Compute the estimated indices γ ∈ R n on Y and the whole covariate matrix X using SIR-QZ.

Let a = 1.

Step 2. Select randomly p 0 components of x and build the corresponding matrix X (a) .

Step 3. Compute the SIR-QZ indices γ(a) ∈ R n based on Y and X (a) .

Step 4. Calculate the linear correlation between the indices γ and γ(a) . Let us denote by ĉ(a) the square of this correlation.

Let a = a + 1.

Repeat N 0 times steps 2-4.

Step 5. Consider the submodels corresponding to the N 1 largest correlations ĉ(a) .

Either the user set ζ ∈]0, 1[ and then gets N 1 = ζN 0 , or the user chose a value for ρ and then N 1 is the number of submodels such that ĉ(a) > ρ.

Step 6. Count the number of occurrences of each component of x in these N 1 submodels. The components that affect y are the ones that have the greater number of occurrences.

For example, in our simulation study, we set N 0 = 10 4 and ζ = 10% to determine the closest N 1 = ζN 0 submodels. In our real data application, we use N 0 = 9 × 10 5 and ρ = 0.9 to select the top N 1 submodels.

Note that choosing p 0 < n allows us to use classical SIR instead of SIR-QZ in Step 3 which significantly improves the computational time. In addition, any SIR approach that provides estimates of the indices (when n ≤ p and n > p) could be used in the whole algorithm instead of SIR-QZ.

A simulation study

In Sections 4.4.2-4.4.3, we presented 4 methods to estimate EDR directions (or indices) and 3 procedures to select which components of x have effects on y. In Section 4.4.4.1, we illustrate them on a single simulated data set. To compare their numerical performances, we then study them on several replications in Section 4.4.4.2. where x and β are p-dimensional vectors defined hereafter and ε ∼ N (0, σ 2 ) with σ = 10 -3 . Let p = 200 and p ⋆ = 20. We choose β = 1 10 (I(1 ≤ p ⋆ ), . . . , I(p ≤ p ⋆ )) ′ , so that p ⋆ is the number of non-null components of β, that is the number of components of x that affect y. We construct x = (x 1 , . . . , x p ) ′ as follows: for j = 1, . . . , p ⋆ , generate σ 2 j from the law U ([0.05, 0.1]) and x j from the law N (0, σ 2 j ). For j = p ⋆ + 1, . . . , p, set

σ 2 j = 12 -⌊(j -1)/p ⋆ ⌋ ⌊(j -1)/p ⋆ ⌋ 2 σ 2 (j-1) mod p ⋆ +1 ,
when mod denotes the modulo operation. Generate then xj from the law N (0, σ 2 j ) and set x j = x (j-1) mod p ⋆ +1 + xj . This ensures that cor(x j , x (j-1) mod p ⋆ +1 ) = ⌊(j -1)/p ⋆ ⌋/12.

Estimation of EDR indices. We simulate an independent and identically distributed sample (X,Y ) of size n = 100 from model (4.29). We plot Y versus the true indices X ′ β in Figure 4.13. We analyze X and Y with the various methods presented in Section 4.4.2.2-4.4.2.4.

• For the RSIR method we evaluate L(s) for s ∈ {0, 10 -10 , 10 -9 , . . . , 10 5 } with 50 bootstrap samples and H = 10. In Figure 4.14(a) we plot the values of Y against the indices provided by the RSIR method. The structure of Figure 4.13 can not be discerned in Figure 4.14(a). The regularization parameter that RSIR provides is equal to 10 5 and thus the RSIR procedure is equivalent to an eigendecomposition of M .

• Concerning the SR-SIR method, the regularization parameter s opt is chosen among values in {10 -10 , 10 -9 , . . . , 10 5 }. The number of iterations of the alternating least square algorithm of SR-SIR is set to 50 and we take H = 10. For this example, we find s opt = 10 3 . In projector on the subspace of R n spanned by X ′ β. More precisely, we have

P = X ′ β(β ′ X X ′ β) -1 β ′ X. (4.30)
Similarly, define P RSIR , P SR-SIR , and P SIR-MP by respectively replacing β by βRSIR , βSR-SIR , and βSIR-MP in (4.30). Let us also define P SIR-QZ by

P SIR-QZ = Īn γ(γ ′ Īn γ) -1 γ′ Īn ,
where Īn is a matrix that centers γ (see Appendix A.2.3). Note that for any a = 0, we also have P SIR-QZ = Īn (aγ)((γ ′ a) Īn (aγ)) -1 (γ ′ a) Īn , which is coherent since if (γ, δ) is a solution of (4.27), then (aγ, 1 a δ) is another solution of this very equation. For a given method m, we use the squared trace correlation between the subspaces spanned by X ′ β and by X ′ βm as a measure of the closeness between these subspaces. It is defined by

R(m) = 1 K Trace(P P m ). (4.31) Notice that if K = 1, R(m)
is the squared cosine of the angle between the vectors X ′ β and X ′ βm . This quality measure belongs to [0, 1] and the higher its value is, the better the indices are estimated.

In Table 4.1, we present values of R(m) for the four considered methods. We also present, in this table, the computational time needed for our Intel Core 2 Quad Q9505 processor to execute these methods. For SIR-QZ, R(m) is clearly higher than for RSIR, SR-SIR and SIR-MP (0.74 versus less than 0.1). Notice also that the computational time is a lot greater for RSIR and SR-SIR than for SIR-MP and SIR-QZ.

Selection of components of x. We rely on the true positive rate (TPR) and on the false positive rate (FPR) to evaluate procedures that find which elements of β are equal to 0. The TPR is the number of selected components of x that actually affect y divided by the total number of components of x that affect y. The FPR is the number of selected components of x that do not affect y divided by the total number of components of x that do not affect y.

For the RSIR selection method the returned p-values are ordered and the components that correspond to the first p-values are selected. For the CSS procedure, the components related to the greatest number of occurrence are selected. We evaluate the results of both methods by selecting the best 10 (resp. 20, 50, 100, 150) components. For the SR-SIR method, we use the different criteria (AIC, BIC and RIC) proposed in [START_REF] Li | Sliced inverse regression with regularizations[END_REF] to determine the appropriate Lasso parameter.

The number of bootstrap samples generated for the RSIR method is set to 10 3 . It is great enough so that increasing it does not improve significantly the quality criterion. This method relies on a regularization parameter which can also be provided by the Algorithm A.1 to find βQZ (see Appendix A.2.2). The estimate βQZ can thus be plugged in the RSIR selection method. For the CSS method, we choose N 0 = 10 4 , ζ = 10% and p 0 = 50. While increasing N 0 could lead to better quality criteria, the computational time is sufficiently large not to choose it greater. The tested values of the Lasso parameter for the SR-SIR method are in {1, 2, . . . , 100}. This algorithm needs an estimate of β in input. We use βSR-SIR but, because of the poor results of Table 4.1 for this estimate, we also take βQZ with H = 10, and the true EDR direction β.

Results of the corresponding TPR and FPR are displayed in Table 4.2. The SR-SIR method performs very well if an accurate estimate of β is provided, while the results are really bad otherwise because no selected component of x has any effect of y. Concerning RSIR and the CSS method, their FPR are similar, but the TPR are greater for the latter. The results for RSIR with β QZ are slightly better than with the full RSIR procedure. The CSS method also seems to need good estimates of the indices since working with the true ones produces better TPR than using γ.

To get more insights about the numerical performances of the various procedures tested in this example, we run them in several replications in the following section.

General behaviors of the estimates over several replications.

We generate 100 samples of size n = 100 from the model (4.29). For each of them, we launch the RSIR, SR-SIR, SIR-QZ and SIR-MP procedures with the previously described parameters, and compute the quality criterion given in (4.31). We display the values of the criterion in Figure 4.15. The trend that is exhibited in Table 4.1 is confirmed in this graphic since the values of R(SIR-QZ) are clearly greater than the others. • RSIR with βRSIR , 1000 bootstrap samples and the following significance levels of the corresponding test: 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,

• SR-SIR with βSR-SIR and τ ∈ {1, 2, . . . , 100},

• The CSS procedure with γ from SIR-QZ and parameters N 0 = 10 4 , ζ = 10% and p 0 = 50 for various number of selected components: 10, 20, 30, 40, 60, 80, 100, 120, 140, 160. The mean ROC curves over the 100 replications are displayed in Figure 4.16. The CSS method outperforms RSIR while SR-SIR provides poor results. Notice that for RSIR, the values of the FPR are close to the chosen levels of test, in spite of the fact that this test is asymptotic.

Remark 15. The simulated model (4.29) may appear complicated. We explain here why we chose it by considering the following one:

y = sign(x ′ β 1 ) log(|x ′ β 2 + 5|) + ε, (4.32)
where p = 200 and ε ∼ N (0, 0.09). In addition, β = (β 1 , β 2 ) satisfies

β 1 = (1, 1, 1, 1,

SIR-QZ and CSS results

We first run SIR-QZ for {H 1 , . . . , H N H } = {2, . . . , 6}. In Figure 4.18(a), we plot the dependent variable Hopx versus the index based on the whole set of SNPs (p = 770). This graphic clearly exhibits a link between the phenotype and the index estimated by a smooth kernel method. In this illustration, this link is almost linear. In contrast, the plot (not given here) of Hopx gene versus the second EDR index does not show any structure. From this graphical diagnostic, it appears that only one EDR direction provides relevant information to explain the variability of the gene expression. Then, in order to find a parsimonious index, we run our CSS selection procedure with N 0 = 900.000. The examined values of p 0 are in {10, 20} while ρ takes value in {0.75, 0.80, 0.85, 0.90}. In Table 4.3, we present the number of selected SNPs for each combination of this two parameters. The threshold used for the selection will be detailed below. We can observe that, not surprisingly, the numbers of selected SNPs increases with p 0 and decreases with ρ. Moreover for a given value of p 0 , we specify, in this table, the number of selected SNPs in common with those selected with ρ = 0.9 (corresponding to the parsimonious model). We also indicate, for a given ρ, the number of SNPs in common with those selected when p 0 = 10 and when p 0 = 20. This table highlights an overlap of 10 SNPs among all the sets of the selected SNPs for the various couple (p 0 ,ρ). Note that, the smallest set contains 11 SNPs when p 0 = 10 and ρ = 0.9 which comforts us about the stability of the CSS procedure.

In eQTL study, it is known that only a few number of SNPs can explain the variation of the gene expression. Thus, from the expertise of the biologists, we decide to select the Figure 4.20 presents the Venn diagram of the sets of SNPs selected by the different approaches. Two SNPs (D14Mit3 and D2Cebr204s17) are selected by the four methods. D14Mit3 (chromosome 14) is clearly the first SNP in the list of the SNPs selected by the CSS procedure (see Figure 4.19) and D2Cebr204s17 (chromosome 2) is at the third position. Moreover, our proposed approach reveals 4 SNPs (D18Ucsf1, Fgg, D2Utr9, D16Rat75) not selected by the other methods. Although three of them are very close to our proposed threshold (red line in Figure 4.19), while the SNP D18Ucsf1 (chromosome 18) is clearly selected by our procedure.

The main advantage of our approach is the opportunity to reveal a non-linear link between the gene expression and a parsimonious index while the other compared approaches are based on a linear model. However, theses methods could treat multiple correlated phenotypes (multiple continuous responses). For example, ESS++ has been used to study the joint variability of gene expression in seven tissues (adrenal, gland, aorta, fat, heart, kidney, liver, skeletal) from inbred line rats [START_REF] Bottolo | ESS ++ : a C ++ objectedoriented algorithm for bayesian stochastic search model exploration[END_REF]). An extension of SIR-QZ for multivariate response is under investigation.

Concluding remarks

Although regularizing the estimated covariance matrix or constraining the optimization problem are natural ways to extend SIR to underdetermined cases (n < p), it may not be clear which one should be chosen and how to set the related parameter for each procedure. For the RSIR method, we illustrated that in such a context, the corresponding parameter should rather be determined with respect to the stability of the linear algebra algorithm (as in SIR-QZ) than with a statistical criterion. Moreover, the SIR-QZ approach introduced in this paper produces better results in simulation than the SR-SIR method that constraint the underlying optimization problem. In addition, the poor performances of SIR-MP suggest that adapting properties of the pseudo-inverse from the functional SIR to our high dimensional context is not well-adapted. We assumed that the dimension K of the EDR subspace was known in our simulation study. While, in the application, an empirical argument was given to determine this dimension, its estimation remains to be done with care.

We also proposed the CSS method that searches which submodels are the most informative to select relevant components of x. This procedure relies on a quality measure of the estimated indices, and it outperforms, in simulation, RSIR and SR-SIR selection procedures which are both based on estimates of β. We thus explain these results by pointing out that, when n < p, only the indices can be properly estimated, in general cases. Note that the space of the submodels may not be browsed optimally by our proposed CSS algorithm. Improvements should be made using optimization techniques such as genetic algorithms.

All the developed methods have been implemented in R language and are available upon request from the corresponding author. An R package is currently under development.

The illustration of genetic data highlights the opportunity of SIR-QZ combined by the CSS procedure to reveal a few number of SNPs which can explain the variability of the expression of the Hopx gene. Then, a non linear link between the gene expression and the parsimonious index could be estimated. The choice of the number of SNPs to keep is still a topic of concern since alternatives to our threshold could be considered.

In genetic datasets, the response variable is often multivariate. For instance, it could represent several phenotypes as in eQTL studies. Some approaches already handle such datasets. Since univariate results for SIR-QZ and the CSS selection procedure are consistent with the other methods presented in our application, it thus appears interesting to extend them to the multivariate case.

Chapter 5 Conclusion

In this document, we exhibit various stochastic models adjusted for some biological issues. A particular attention was given to include knowledge and needs from the tackled problems into the statistical analyses. When several methods have been developed for a similar goal, another concern was to compare them using theoretical arguments or simulation studies. The procedures we consider were made accessible by implementing them in a user-friendly way.

The model of Section 3.1 was designed to handle reasonable assumptions about the valvometric signal. Various models exist to detect states in a stochastic process defined as a piecewise function on which a noise is added. For valvometric signals, we would rather assume that the underlying function is continuous because the distance between both parts of a shell of an animal does not vary with an infinite speed. In addition, two consecutive points in the signal do not seem to be independent even with respect to the current state of the hidden process. For example, if a point belongs to a microclosing, the following one may not follows the same distribution as if it does not. This can be observed in Figure 3.16. The procedure we present to estimate the hidden process does not thus need independent observations but only equally spaced points.

Statistical needs related to genetic datasets, such as multiple testing, are well-known and widely studied. That is why we worked in an existing framework called FAMT which can only handle linear links. Because we anticipate that non-linear ones may also be of interest, we gave in Section 4.3 a method to apply SIR in such a multivariate context. Due to the cost of DNA sequencing, the sample size of a transcriptomic dataset may be quite small. That is why we also provided some modifications of SIR adapted to cases when n < p, in Section 4.4.

Finding a threshold to separate the open state from the closed one in the valvometric data requires to estimate a probability density function f with a single anti-mode located in θ. Various alternatives are available for this purpose. We showed in the simulation study of Section 3.2.3 that computing Polonik's estimator may not produce good enough estimates of θ. When f have bounded support, the kernel-based estimators with Gaussian kernel seems to succeed in this task whether h SJ or h crit is chosen. Appendix A.1 gives a possible reason why using h crit with the uniform kernel may provide inconsistent estimates of θ.

Multivariate sliced inverse regressions can be performed with various procedures like pooled marginal slicing or k-means inverse regression. The results of the simulation study of Section 4.3.5 let us think that these methods manage to properly estimate EDR spaces, as so do MSIR and wMSIR. The latter ones can however straightforwardly be paired with a procedure to detect clusters of components in the multivariate response variable.

When the sample size is smaller than the number of explanatory variables and in the particular context of model (4.29), SIR-QZ outperforms RSIR, SR-SIR and SIR-MP when trying to estimate a basis of the indices X ′ β. The CSS procedure, together with SIR-QZ, also achieves a better selection of the explanatory variables than the full RSIR and SR-SIR methods, for this model.

A lot of theoretical results in statistics are given when n → ∞. Statistics is also a privileged science when considering high-dimensional variables. It thus makes sense to have to handle large data sets to apply statistical methods. Computer programs which help to do so are then welcome. SQL databases are appropriate structures for storing and easily accessing such data sets. Free programs exist to manage these databases. We then imported most real data we can into databases in the analyses presented in this thesis. We also explained how to retrieve data from the SQL databases with the R software in Section 3.4.

Executing R functions built by statisticians and sending SQL requests to another server may take a lot of time for occasional R users. Deploying a graphical user interface (GUI) with an R package may circumvent these difficulties. We created a GUI to estimate the probability density f from valvometric data (see Section 3.4.3) and another one in order to apply functions from the FAMT package on the IMMORTEEL transcriptomic data (Section 4.2.3).

Concerning the estimator f n described in equation (3.6) with h n = h crit , further work would be interesting. The asymptotic behaviour of the bandwidth h crit should be studied in order to know whether or not it goes toward 0 and if nh 2 crit → ∞, when n → ∞. This is important in order to have Assumptions 5 and then convergence of the position of the unique local minimum of f n toward θ.

For a kernel density estimator f n of f , with a random bandwidth, [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF] and [START_REF] Devroye | A Course in Density Estimation[END_REF] gave the convergences of f nf ∞ and of f nf 1 . For a kernel density estimator f n with a real bandwidth, there exist results about f nf 2 . To our knowledge this is not the case for random bandwidths. For these bandwiths, proving that f nf 2 converges toward 0 or explaining why it is not possible may be an exciting challenge. It could for example bring new arguments in the comparison between Polonik's estimator and a kernel estimator relying on the bandwidth h crit .

In the example of Section 3.2.3, we actually found that for a kernel density estimator with the Gaussian kernel, choosing the bandwidth h crit,k with a large integer k may produce smaller IAE values than taking the bandwidth h crit = h crit,N (f ) . It also creates spurious modes in fK,h crit,k . In his test of the existence of modes, [START_REF] Minnotte | Nonparametric testing of the existence of modes[END_REF] proposed to replace some modes with constant parts. Using this smoothing technique, we could have for all k ≥ N (f ) an estimator based on fK,h crit,k with only N (f ) modes. Theoretical properties of this estimator remain to be found.

In Section 3.1 we classified some oysters into two groups. This procedure only relies on their transitions between their open and their closed state. Their microclosings may however furnish additional information about their health status. Considering their amplitude and the time between two consecutive microclosings, we could for example study the bivariate distribution of these variables like in [START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF]. This requires a microclosing detection technique, which could be built from Section 3.3.3.

The R package we created to access valvometric data sometimes need a significant amount of computational time. To improve it, we could create another index in the SQL database which would consider the location of the studied animal, the date and the bivalve's identification number, in this order. We also could implement a cache for the valvometric signal. For example, when a user studies data from a given day, the program could silently download signals related to the next or to the previous ones in order to display them faster if the user ask for them.

The FAMT model (4.5) is a particular case of model (4.9), with K = 1. We can then apply the MSIR method on the IMMORTEEL dataset in order to find clusters in the set of selected genes. If for j 1 ∈ {1, . . . , q} and i 1 ∈ {2, . . . , p}, we have equation (4.6) with j = j 1 and i = i 1 and if we have it too for j = j 2 and i = i 2 where j 2 ∈ {1, . . . , q} and i 2 ∈ {2, . . . , p}\{i 1 }, then y (j 1 ) and y (j 2 ) should not appear in the same cluster.

In model (4.5), let Y := y (1) , . . . , y (q) , and β := β (1) , . . . , β (q) . Let also ε := ε (1) , . . . , ε (q) and B := b (1) , . . . , b (q) . If the rank of β is equal to p, because p < q for the IMMORTEEL dataset, we have

X ′ = Yβ + -Z ′ Bβ + -εβ + ,
where β + is the Moore-Penrose pseudoinverse of β. Let β + =: β + 1 , . . . , β + p . We can write for all i ∈ {1, . . . , p}

X i = Yβ + i -Z ′ Bβ + i -εβ + i .
We have n = 44 and q = 17, 866 which are appropriate values to execute SIR-QZ and the CSS procedure on the response variable X i and on the matrix of explanatory variables Y.

While SIR was introduced in 1991, new extensions can still be proposed nowadays. The underlying model (4.18) could for example be modified so that for all k ∈ {1, . . . , K}, β k ∈ N p . This might make the null elements of β k easier to detect, when trying to select relevant explanatory variables. This approach was studied in [START_REF] Sabatier | Extensions of simple component analysis and simple linear discriminant analysis using genetic algorithms[END_REF] and [START_REF] Vines | Simple principal components[END_REF] for the principal component analysis.

Finally, it seems to us that the statistical models studied in this document may also be useful for other biological problems and in other fields than biology. Continuous variations in a signal in which a process with a finite number of states is hidden may not only happen for valvometric measures. Graphics in [START_REF] Watwood | Deep-diving foraging behaviour of sperm whales (physeter macrocephalus)[END_REF] suggest that cetaceans' dive could also be studied with the framework of Section 3.1. SIR and FAMT are suitable methods when aiming at finding a few useful explanatory variables among a lot of considered ones. In marketing, we may find the reasons why people buy a given product with these statistical tools. Theorem 15. For any probability density function f of X, let fK,h crit,k be the estimator of f when K is the uniform kernel with h crit,k given in (A.1). Then we have h crit,k increasing with n, for all k ∈ N.

The proof is given in Appendix A.3.9.

To illustrate this theorem, we consider again the sample (3, 6.5, 6, 1.5, 5) and every subsamples made of the first n elements of (3, 6.5, 6, 1.5, 5) for n ∈ {1, . . . , 5}. For each subsample, we compute the function h → N ( fK,h ), where K is the uniform kernel. In Figure A.2, we display ranges of values of h and n for which N ( fK,h ) is equal to a given number. The increase of N ( fK,h ) with n can be observed for small values of h. In the general case, this feature implies that for any k ∈ N, h crit,k also increases with n.

Theorem 15 means that we can not have that h crit,k goes toward 0 in probability for the uniform kernel. Thus, the proof of the pointwise convergence of fK,h crit,k toward f , as given in [START_REF] Futschik | On the consistency of kernel density estimates under modality constraints[END_REF], does not hold for this kernel. It is also impossible to use the work of [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF] and [START_REF] Devroye | A Course in Density Estimation[END_REF] to show the corresponding L 1convergence and uniform convergence. In addition, the simulation study from Section 3.2 shows that fK,h crit,k is not an accurate estimator of f , when K is the uniform kernel. For these reasons, we recommend not to use fK,h crit,k with this kernel in practice. (3, 6.5, 6, 1.5, 5). Proofs of ΣBB ′ Σ and Σ associated with K non-null eigenvalues. Thus, V contains the K eigenvectors of Σ -1 (ΣBB ′ Σ) = BB ′ Σ which are associated with K non-null eigenvalues.

In addition, we have that Span(ΣBB ′ Σ) = Span(ΣB)

because Span(ΣBB ′ Σ) ⊂ Span(ΣB) and dim(Span(ΣBB ′ Σ)) = K. Since Σ is invertible, this implies that Span(BB ′ Σ) = Span(B).

Finally, we have under model (4.9) and assumption (LC) that for all j ∈ {1, . . . , q}, Span(B) = Span(B (j) ) using univariate SIR theory. We then have Span(B) = Span(B) and then Span(V) = Span(B).

A.3.3 Proof of Theorem 10

For each component y (j) of y and under the assumptions (LC), (A1)-(A3), from univariate SIR theory of [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], each estimated EDR space B (j) converges to B (j) at root n rate: that is, for j = 1, . . . , q, B (j) = B (j) + O p (n -1/2 ). It follows that B = B + O p (n -1/2 ). Since Σ = Σ + O p (n -1/2 ), we get B B ′ Σ = BB ′ Σ + O p (n -1/2 ). Therefore, the eigenvectors associated with the largest K eigenvalues of B B ′ Σ converge to the corresponding ones of BB ′ Σ at the same rate: vk = v k +O p (n -1/2 ) for k = 1, . . . , K. Consequently, the estimated EDR space Span( V) converges to Span(V) at rate √ n. Since Span(V) = Span(B), the estimated EDR space converges to the true one in probability.

A.3.4 Proof of Theorem 11

We can proceed analogously to the proof of Theorem 10. It is sufficient to show that W = W + O p (n -1/2 ). Since, for each component y (j) of y and under assumptions (LC) and (A1)-(A3), from univariate SIR theory of [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF], each estimated eigenvalue λ(j)

k converges to λ (j)
k at rate √ n, we have π(j) = π (j) + O p (n -1/2 ) and π⋆ = π ⋆ + O p (n -1/2 ).

We thus get W (j) = W (j) + O p (n -1/2 ) and W = W + O p (n -1/2 ). Consequently, B W B ′ Σ = BWB ′ Σ + O p (n -1/2 ) and the eigenvectors associated with the largest K eigenvalues of B W B ′ Σ converge to the corresponding ones of BWB ′ Σ at the same rate: ṽk = ṽk + O p (n -1/2 ) for k = 1, . . . , K. Therefore, the estimated EDR space Span( V) converges to Span( V) at root n rate. Since Span( V) = Span(B), the estimated EDR space converges to the true one in probability.

A.3.5 Proof of Lemma 12

For j = 1, . . . , q, recall that P B(j) , Σ is the Σ-orthogonal projector onto Span( B(j) ). According to Theorem 4 of [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF], we have that P B(j) , Σ = P B 1 ,Σ + O P (n -1/2 ), for A.3.9 Proof of Theorem 15

First, note that when K is the uniform kernel, for some h > 0, we can find N ( fK,h ) by counting the number of variations of sign of the following function where ε is chosen in a way that ensures that ∀(i, j) ∈ {1, . . . , card(A h )} × {1, . . . , card(B h )} , (a h,(i)b h,(j) ) ∈] -∞, 0]∪]ε, ∞[, in order to obtain a unique value of g ′ h,ε (x) for each x. The aim of g ′ h,ε is to mimic the derivative of fK,h . It seems easier to use than dirac functions involved in f ′ K,h . Besides, one can see that N (g h,ε ) = N ( fK,h ), using the fact that Proposition 14 is valid for g h,ε . That is why the number of variations of sign of g ′ h,ε is equal to 2N ( fK,h ) -1. Let C ε,n = c h,ε,(i) i∈{1,...,w} be the ordered sequence made of the sets a h,(i) -ε 2 i∈{1,...,card(A h )} and b h,(i) + ε 2 i∈{1,...,card(B h )} .

Let d h,ε,(i) = 1 c h,ε,(i) ∈ a h,(i) -ε 2 i∈{1,...,card(A h )}

-1 c h,ε,(i) ∈ b h,(i) + ε 2 i∈{1,...,card(B h )} , and D ε,n = d h,ε,(i) i∈{1,...,w} . Every interval where g ′ h,ε (x) = 0 is represented by a c h,ε,(i) , then the number of variations of sign is the same for g ′ h,ε and for D ε,n . We write v(D ε,n ) for the number of variations of sign of D ε,n like [START_REF] Schoenberg | On pólya frequency functions II: Variation diminishing integral operators of the convolution type[END_REF] did in his article. Now, we prove that v(D ε,n ) ≥ v(D ε,n-1 ), for n > 1. This property is satisfied if D ε,n-1 = JD ε,n where J is a totally positive matrix, following [START_REF] Schoenberg | On pólya frequency functions II: Variation diminishing integral operators of the convolution type[END_REF]. To define J, we first focus on the case where the last point in the sample is different from the others. This means that if Ω is our sample space, we define Ω 1 as: Ω 1 = {ω : ∀i ∈ {1, . . . , n -1} , X i (ω) = X n (ω)} .

We remark that, when our sample comes from ω ∈ Ω 1 , D ε,n-1 is constructed by removing two points in D ε,n . These points correspond to c h,ε,(γ 1 ) = X nh -ε 2 and c h,ε,(γ 2 ) = X n + h + ε 2 . This is why we have :

J =                 I γ 1 -1 0 . . . 0 • • • • • • 0 . . . 0 • • • 0 . . . • • • 0 0 • • • . . . 0 • • • 0 . . . 0 I γ 2 -γ 1 -1 0 . . . 0 • • • 0 . . . • • • 0 0 • • • . . . 0 • • • 0 . . . 0 • • • • • • 0 . . . 0 I w-γ 2                

Figure 1 . 1 :

 11 Figure 1.1: Un exemple de données valvométriques, pour une huître de la jetée d'Eyrac, le 27 février 2007

Figure 2 . 1 :

 21 Figure 2.1: An example of valvometric data, for one oyster at Eyrac pier on the 27th of February 2007

Figure 3 .

 3 Figure 3.1: Two electrodes stuck on an oyster. Picture taken by the MolluSCAN Eye team.

Figure 3 .

 3 Figure 3.2: Example of trajectory of the four-state renewal process for modeling oyster's behavior

. 4 )

 4 Proof. Azaïs et al. have shown that the multiplicative intensity model makes sense for the observation of the sequence (Z k , S k ). The proof of (3.4) under the conditions (3.2) and (3.3) can be found in Andersen et al. (1993), Theorem IV.1.1, when the multiplicative intensity assumption is verified. Lemma 2. Let x ∈ E. Assume that λ(x, •) is a locally integrable function. Then, for any t ≥ 0, the conditions (3.2) and (3.3) are verified. Proof. Given in Appendix A.3.1.

Figure 3 . 3 :

 33 Figure 3.3: Fleming-Harrington estimates of the survival function with respect to t for four oysters: (solid line) oyster of the Eyrac pier between 1st and 8th of May 2007; (dashed line) oyster of Locmariaquer between 1st and 8th of October 2011; (dotted line) oyster of Santander between 1st and 8th of December 2011; (mixed line) oyster of Eyrac pier between 12th and 19th of May 2009. 3.3(a) Survival function estimation for x = C; 3.3(b) Survival function estimation for x = O.

Figure 3 .

 3 Figure 3.4: Fleming-Harrington estimates of the probability that the first return time to state T OC is greater than t for four oysters: (solid line) oyster of the Eyrac pier between 1st and 8th of May 2007; (dashed line) oyster of Locmariaquer between 1st and 8th of October 2011; (dotted line) oyster of Santander between 1st and 8th of December 2011; (mixed line) oyster of Eyrac pier between 12th and 19th of May 2009.

Figure 3 . 5 :

 35 Figure 3.5: Projection of classified oysters on the first two principal component axes which display 97% of the total variability. The label of each oyster corresponds to its location: E: Eyrac, L: Locmariaquer, S: Santander, (d): Dying oysters, (t): Oyster affected by quick temperature variations.

  (a) Estimations fK,hcrit with Gaussian kernel (dashed line) and fK,hcrit with uniform kernel (mixed line). (b) Estimations fP (long dashed line) and fK,h SJ (dotted line).

Figure 3 . 6 :

 36 Figure 3.6: Beta mixture density f 1 (solid line), and various estimations; for n = 1600.Vertical lines: positions of the minimum of local minima of each plotted density. In Figure3.6(b), vertical axis is broken between 2.2 and 2.8, between 3.5 and 5 and between 5.5 and 12.5 to be able to clearly see the shape of all curves.

  (a) Kernel density estimator with h crit and the uniform kernel. (b) Kernel density estimator with h crit and the Gaussian kernel. (c) Kernel density estimator with h SJ and the Gaussian kernel. (d) Polonik's estimator.

Figure 3 . 7 :

 37 Figure 3.7: IAE for the beta mixture model (3.12) and various density estimators.

  (a) Estimator related to fK,hcrit with the uniform kernel. (b) Estimator related to fK,hcrit with the Gaussian kernel.(c) Estimator related to fK,h SJ with the Gaussian kernel.(d) Estimator related to fP .

Figure 3 . 8 :

 38 Figure 3.8: Estimations of the position of the local minimum of the density of the beta mixture model.

  (a) Estimations fK,hcrit (dashed line) and fK,h SJ (dotted line). (b) Estimations fK,hcrit,2 (mixed line) and fK,hcrit,7 (long dashed line).

Figure 3 . 9 :

 39 Figure 3.9: Asymmetric claw density f 2 (solid line) and various estimations; for n = 1600 and the Gaussian kernel K. Vertical lines: positions of the minimum of local minima of each plotted density.

  (a) Kernel density estimator with h crit and the Gaussian kernel. (b) Kernel density estimator with h SJ and the Gaussian kernel.

Figure 3 .

 3 Figure 3.10: IAE for the asymmetric claw model (3.13) and various density estimators.

  (a) Estimator related to fK,hcrit with the Gaussian kernel. (b) Estimator related to fK,h SJ with the Gaussian kernel.

Figure 3 .

 3 Figure 3.11: Estimations of the position of the local minimum of the density of the asymmetric claw model.

Figure 3 .

 3 Figure 3.12: Density f 3 of a mixture of t-distributions (solid line). Estimations fK,h crit with Gaussian kernel (dashed line) and fK,h SJ (dotted line), for n = 1600. Vertical lines: positions of the minimum of local minima of each plotted density.

  (a) Kernel density estimator with h crit and the Gaussian kernel. (b) Kernel density estimator with h SJ and the Gaussian kernel. (c) Polonik's estimator.

Figure 3 .

 3 Figure 3.13: IAE for the mixture (3.14) of t-distributions and various density estimators.

  (a) Kernel density estimator with h crit and the Gaussian kernel. (b) Polonik's estimator.

Figure 3 .

 3 Figure 3.14: Estimations of the position of the local minimum of the density of the mixture of t-distributions with various estimators.

  (a) Oyster assumed to be healthy.(b) Oyster that dies the week that follows the day when measurements were taken.

Figure 3 .

 3 Figure 3.15: Density estimations with fK,h crit,2 of opening amplitudes of two oysters. Crosses indicate local extrema of the estimated densities.

Figure 3 .

 3 Figure 3.16: Valvometric measures from an oyster living at Eyrac on May 8, 2007 from 00:34:10 to 00:40:59.

  , . . . , 0, -47.58, 43 0, . . . , 0, -19.60, 33 0, . . . , 0, -205.99, 134.03, 23 0, . . . , 0, 30.63, 11 0, . . . , 0).

Figure 3 .

 3 Figure 3.17: Wavelet coefficients computed from {ỹ k } k∈D J (left) and {ẙ k } k∈D J (right). Dashed line: threshold δ h = σ 2 log(n log(n)/ log(2))/ √ n.

Figure 3 .

 3 Figure 3.18: Valvometric signals reconstructed after hard thresholding of the wavelet coefficients of U j 0 computed from {ỹ k } k∈D J (left) and {ẙ k } k∈D J (right), using j 0 = ⌊log(log(n))/ log(2)⌋ + 1 and δ h = σ 2 log(n log(n)/ log(2))/ √ n.

  For another part of this signal during which the oyster is open, let b := max k∈D J {ỹ k } and c := b-max k∈D J {ỹ k : ỹk = b}.

Figure 3 .

 3 Figure 3.19: Microclosing detection using the first 32, 769 points of Figure 2.1. Crosses are locations where they are detected.

  The data.frame objects are tables and can contain both strings and real numbers, among other kind of variables. The read.table function creates such an object from some data saved in a text file. Let my.d.f be a data.frame with two lines and three columns. The command dimnames(my.d.f)[[1]] <-c("Fst.line", "Snd.line") changes the names of the lines while dimnames(my.d.f)[[2]] <-c("Fst.col.", "Snd.col.", "Trd.col.") does the same for the columns. Then, my.d.f[1,] and my.d.f["Fst.line",] both return the first line of my.d.f, while the commands in Algorithm 3.1 display its second column.

>

  my.d.f[,2] > my.d.f[,"Snd.col."] > my.d.f$Snd.col. > attach(my.d.f) > Snd.col. Algorithm 3.1: Different methods to obtain the second column of a data.frame > my.function.double <-function(x) { x^2 } > my.function.character <-function(a) { rep(a,2) } > my.function <-function(object){UseMethod("my.function")} > my.function(3) [1] 9 > my.function("test") [1] "test" "test" Algorithm 3.2: An example of object-oriented programming using S3 the classes of objects for which the method my.function exists. It returns an explicit message otherwise.

  Figure 3.20: Samples of pairs (X, Y ) from model (3.21)

  mu)^3 + rnorm(n * nMu, 0, 1.5) > par(mar=c(5,5,1,1)) > plot(X,Y,cex.axis=1.5, cex.lab=1.5) <Figure 3.20(a)> > trellis.par.set("strip.shingle", + list(col=rgb(0.5,0.5,0.5))) > trellis.par.set("strip.background",list(col=0)) > xyplot(Y~X|mu, col=1) <Figure 3.20(b)> Algorithm 3.3: Scatter plots with and without the lattice package $ mysql -h 0.0.0.0 -u user.name -p Algorithm 3.4: Reaching a MySQL server with mysql mysql> CREATE DATABASE myDb; Algorithm 3.5: How to create an SQL database mysql> CREATE

Figure 4

 4 Figure 4.1: Explanatory variables projected onto the first two components of the PCA of X

Figure 4

 4 Figure 4.2: A graphical user interface to choose the explanatory variables to take into account

Figure 4

 4 Figure 4.4: A selection is performed on genes using the thresholds γ and γ ⋆

  Figure 4.5: Comparison of estimators of B on 100 samples from model (4.15) with n = 100 and p = 20. The line in Figure 4.5(b) corresponds to the first bisecting line.

  Figure 4.6: Study of the contribution of y (3) to V for 100 samples generated from model (4.15) with n = 100 and p = 20.

Figure 4

 4 Figure 4.7: Boxplots of quality measure for the wMSIR method, for 100 samples generated from model (4.15) (Figure 4.7(a)) or from model (4.16) (Figure 4.7(b)) with various values of n and p.

Figure 4 . 9 :

 49 Figure4.9: Values of rj,j ⋆ in various shades of grays for (j, j ⋆ ) ∈ {2, . . . , q}×{1, . . . , j -1} and for a sample of size n = 1000 generated from model (4.17) with θ 1 = θ 2 = 1 and p = 20.

  Figure 4.10: Clustering analysis for a sample of size n = 1000 generated from model (4.17) with θ 1 = 1, θ 2 = 0 and p = 20.

  Figure 4.11: Values of rj,j ⋆ in various shades of grays for hyperspectral data.

Figure 4

 4 Figure 4.12: Clustering of the components of y and estimates ṽ1 = [ ṽ1,1 , ṽ1,2 , ṽ1,3 ] ′ of the EDR direction for each cluster.

  Figure 4.13: Plot of Y versus X ′ β generated from model (4.29), with n = 100 and p = 200. The horizontal scale was standardized.

  Figure 4.14(b), we draw the values of Y against the estimated indices X ′ βSR-SIR . The points of this graphic do not form the same shape as the points in Figure 4.13. m RSIR SR-SIR SIR-QZ SIR-MP R(m) Quality measure and computational time of various estimates of the indices X ′ β for the simulated sample of size n = 100, with p = 200.

Figure 4 .

 4 Figure 4.15: Boxplots of 100 values of R(m) for various estimation procedures m and samples of size n = 100 generated from model (4.29) with p = 200.

  . . . , 0, 1, 1, 1, 1, p-20 0, . . . , 0) ′ .

Figure 4 .

 4 Figure 4.17: Boxplots of 100 values of R(m) for various estimation procedures m and samples of size n = 50 generated from model (4.32) with p = 200..

Figure 4 Figure 4

 44 Figure 4.18: Plots of the dependent variable Hopx versus the index based on the whole SNPs (Figure 4.18(a)) and on the 10 selected SNPs (Figure 4.18(b)). The linear correlation between these two indices (evaluated on n = 29 rats) is 0.956. The dotted (red) line is the estimated linear model, the solid (blue) line is the kernel estimate of the link function with a bandwidth chosen by cross-validation.

Figure 4

 4 Figure 4.20: Venn diagram of the sets of SNPs selected by Lasso, sPLS, ESS++ and SIR-QZ (combined with CSS) approaches

Figure A. 1 :

 1 Figure A.1: Estimation of a density with the estimator fK,2 and the uniform kernel K (solid line). Solid circles are special points of this estimated density while brackets indicate a difference between the limit and the value of fK,2 . Crosses represent the underlying sample.

Figure A. 2 :

 2 Figure A.2: Evolution of N ( fK,h ) with respect to the bandwidth h and to the sample size n, for the uniform kernel K. The considered samples are the first n values of(3, 6.5, 6, 1.5, 5).

  x ∈ [a h,(i)ε, a h,(i) [, ∀i ∈ {1, . . . , card(A h )} , -1 for x ∈ [b h,(i) , b h,(i) + ε[, ∀i ∈ {1, . . . , card(B h )} , 0 elsewhere,

  = x for only one point for each interval [T k , T k+1 ], for values of k such that Z k ∈ {T OC , T CO }. This leads to the second assumption. Note that because b and c can be arbitrarily close from each other, Assumptions 3 are mild. Provided theses hypotheses, the mean value theorem enables us to show that,

	.1. The
	first one expresses the fact that we are looking for a threshold that separates out values
	of G(t) when the oyster is open or closed. Note that for x ∈ [a, d], f (x) can be linked to the points t ∈ [0, T max ] such that G(t) = x. To compute easily f (x) when x ∈ [b, c] we assumed that G(t)

  Such actions can be remotely done using the Internet. The datasets are then stored in a MySQL server. A MySQL client can be used to access it. If MySQL is included in a commercial product, it is associated with a proprietary license. Otherwise, it is under GPL (General Public License). PostgreSQL and SQLite are fully free alternatives of MySQL. MySQL can be downloaded from http://dev.mysql.com/downloads/mysql/.From text files to an SQL database. Many tutorials about MySQL are available on-line, we then only focus here on a few important commands. Algorithm 3.4 allows Linux users to connect to a MySQL server from a command-line interface. In this line of code, 0.0.0.0 has to be replaced with the IP address of the MySQL server and user.name should be the login the network administrator in charge of the server gave to the user. A password was probably assigned to the latter and he should enter it after typing Algorithm 3.4. Once someone is connected, he can create a new database called myDb using Algorithm 3.5. SQL databases are basically groups of tables. Algorithm 3.6 creates a table myTable with three columns (intCol, stringCol and realCol) that can respectively contain integers, character strings and real numbers. Thanks to Algorithm 3.7, we can thus transfer data to the database from a text file myTable.txt which contains three columns with the appropriate types of variable. Note that the option --local can be added if MySQL does not have sufficient permissions to handle the file myTable.txt.

  for mercury, cadmium, lead, copper and zinc and in[START_REF] Durrieu | Goldmining and mercury contamination of the piscivorous fish Hoplias aimara in French Guiana (Amazon basin)[END_REF] for mercury. PCBs are toxic molecules which could formerly be found in electric devices such as transformers. DDT, dieldrin and lindane are very dangerous insecticides. Living beings can change the first one into DDE and DDD. These molecules are then called metabolites. They are toxic too. HCB is a fungicide. PBDE-47 was used to build fireproof objects. Like PCBs, DDT, dieldrin and lindane, tetrabromodiphenyl ethers and HCB are listed as persistent organic pollutants by a United Nations treaty. Chlordane is another insecticide which belongs to this list. It is a mixture of several components and trans-nonachlor is one of them.

  Table 4.2: True positive rate (TPR), false positive rate (FPR) and computational time of various methods run on the simulated sample of size n = 100, with p = 200, to determine which components of β of model (4.29) are not null.

				m		RSIR	CSS
				Estimates	βRSIR	βQZ	γ	x ′ β
				Quality criteria TPR FPR TPR FPR TPR FPR TPR FPR
	Number	of selected	components	10 20 50 100 150	0.05 0.05 0.10 0.04 0.10 0.04 0.20 0.03 0.15 0.09 0.25 0.08 0.35 0.07 0.40 0.07 0.30 0.24 0.45 0.23 0.60 0.21 0.70 0.20 0.40 0.51 0.70 0.48 0.75 0.47 0.90 0.46 0.85 0.74 0.85 0.74 0.80 0.74 0.95 0.73
	Computational time (s)	120.61	120.64	328.73	321.04
					m	SR-SIR
				Estimates	βSR-SIR	βQZ	β
				Quality criterion TPR FPR TPR FPR TPR FPR
				Selection with AIC	0.00 0.13 0.00 0.16 0.90 0.00
				Selection with BIC	0.00 0.03 0.00 0.07 0.65 0.00
				Selection with RIC	0.00 0.03 0.00 0.04 0.65 0.00
			Computational time (s)	35.53	40.88	33.29

Remerciements Acknowledgements

This work was supported by the University of Bordeaux, the ASPEET project from the University of South Brittany, the Portonovo project and the Centre National de la Recherche Scientifique (UMR EPOC CNRS 5805, UMR CNRS 6205 and UMR CNRS 5251). We are grateful to Jean-Charles Massabuau, Pierre Ciret and Damien Tran from the EPOC laboratory, for the valvometric data and helpful discussions.

Acknowledgements

The authors would like to thank Jean-Charles Massabuau and Pierre Ciret from the EPOC laboratory (University of Bordeaux, UMR CNRS 5805) for the valvometric data they provide us, and for helpful discussions. They are also grateful to Kristin Clements for improvements on the language of this manuscript. Computer time for this study was partially provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. The authors also acknowledge the anonymous referee for useful comments that lead to several ameliorations of this article.

Acknowledgments

The authors are grateful to Efstathia Bura for the Minneapolis schools data she provided us with. They also thank the associated editor and both anonymous referees for their useful comments that lead to several ameliorations of this article.

Acknowledgements

The authors wish to gratefully thank Marie Chavent for the first discussions that led to this work. They also acknowledge the editor and both anonymous referees for their helpful comments which lead to significant enhancements of this article.

A new multivariate SIR approach

The population version of the proposed MSIR approach is first described in Section 4. 3.3.1. Then, the corresponding sample version is given in Section 4. 3.3.2 and some asymptotic results are derived in Section 4. 3.3.3. We then modify MSIR to handle a weighting in the components of y in Section 4. 3.3.4. Methods to choose K and H are discussed in Section 4. 3.3.5. Finally, procedures to withdraw or cluster components of y are detailled in Section 4.3.4.

Population version

Let us assume that the dimension K of the EDR space is known. Let P M,Σ be the Σorthogonal projector on the linear subspace spanned by the columns of a p×K matrix M. A proximity measure between two projectors P M 1 ,Σ 1 and P M 2 ,Σ 2 is given by the squared trace correlation:

for full column rank matrices M 1 and M 2 . For j = 1, . . . , q, let B (j) be a p × K matrix containing a Σ-orthonormal basis of the EDR space from the marginal regression of y (j) given x. Let D be a p × K matrix such that D ′ ΣD = I K , where I K is the identity matrix of order K. Let Q(D, B (1) , . . . , B (q) ) denote the following proximity measure between Span(D) and the q marginal EDR spaces Span(B (1) ), . . . , Span(B (q) ): Q(D, B (1) , . . . , B (q) ) := 1 q q j=1 r D, Σ, B (j) , Σ . (4.12)

This measure takes its values in [0,1]. Note that Span(D) = Span(B (1) ) = • • • = Span(B (q) ) implies Q(D, B (1) , . . . , B (q) ) = 1. The closer to one is this measure, the closer to the q marginal EDR spaces is the linear subspace Span(D).

Let us now consider the following optimization problem: 1) , . . . , B (q) ), (4.13)

where M Σ is the set of Σ-orthogonal p × K matrices. A solution of (4.13) is given by the following theorem.

Theorem 9. Under model (4.9) and assumption (LC), the p × K matrix V is formed by the eigenvectors v 1 , . . . , v K associated with the K non-null eigenvalues of BB ′ Σ where B is the p × (Kq) matrix defined as B := [B (1) , . . . , B (q) ]. Moreover we have Span(V) = Span(B) where Span(B) is the EDR space.

The proof is given in Appendix A.3.2. From Theorem 9, we can estimate a basis of the EDR space based on estimators of matrices B and Σ. This is the goal of the next subsection. 11) = ε (11) , y (12) = ε (12) , (4.17)

based on model (4.15), where (θ 1 , θ 2 ) ∈ [0, 1] and for j = 1, . . . , 12, ε (j) ∼ N 1 (0, 1).

A model with 6 clusters. We first choose θ 1 = θ 2 = 1 which produces the l = 6 following clusters: 2) , y (3) , y (4) , y (5) , y (6) , y (7) , y (8) , y (9) , y (10) , y (11) , y (12) .

We plot, in Figure 4.9, values of rj,j ⋆ for (j, j ⋆ ) ∈ {2, . . . , q} × {1, . . . , j -1}, computed from a sample of size n = 1000. Darker squares correspond to values of rj,j ⋆ close to 1. Thus, Figure 4.9 leads to cluster components y (1) and y (2) together as well as components y (4) , y (5) and y (6) . Components y (7) and y (8) can also be grouped together. It is tempting to put y (3) with y (1) and y (2) and to group y (9) with y (8) and y (7) because the related squared trace correlations are high but r3,9 is approximately as high as them. Note also that for j ∈ {1, 2, 3}, we have r3j,3j-1 < r3j-2,3j-1 and r3j,3j-2 < r3j-2,3j-1 . Recalling that

Clearly, the rank of the p × p matrix Σ is at most equal to n -1 since n i=1 xi = 0 p where 0 p stands for the null vectors of R p . From this remark on the rank of Σ, this matrix is singular when n < p. Moreover, it is also often ill-conditioned when n ≈ p.

Therefore, SIR, SIR-II or SAVE methods only work well when the sample size n is greater than the dimension p of the covariate x, but naturally fail when n < p. In this underdetermined case, the standard estimate of Σ is not invertible even if the components of x are independent.

In the following, we only focus on the SIR approach. We describe it in Section 4.4.2.1 when n > p. The goal of this paper is then twofold. On one hand, we present methods to tackle the issue n < p. On the other hand, we also provide procedures in order to select which components of x have an effect on y.

In Section 4.4.2.2, we consider two different regularizations added to the SIR method, proposed by [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] and [START_REF] Li | Sliced inverse regression with regularizations[END_REF], to find EDR estimates when n < p. Moreover, the SIR method can be seen as a generalized eigenvalue problem and linear algebra algorithms exist to solve this kind of problem without requiring any matrix inversion. The QZ algorithm (see [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF] for instance) is one of them and will be used in the SIR context in Section 4.4.2.3. In Section 4.4.2.4, we also adapt an approach introduced in functional sliced inverse regression (i.e., when x is an explanatory functional variable), based on the Moore-Penrose pseudo-inverse.

Concerning the selection of useful predictors in the indices, [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] use a chi-square test to find which components of x affect y, while the approach of [START_REF] Li | Sliced inverse regression with regularizations[END_REF] relies on a Lasso penalization (Section 4.4.3.1). In Section 4.4.3.2, we propose another procedure. We choose randomly some submodels (i.e., using a number p 0 < p of components of x) and we measure how close they are from the initial one with all the p components of x. The latter model is thus taken as a benchmark. Components of x that appear the most in submodels that are the closest to the benchmark are kept. We naturally consider that the other components of x do not affect y.

In Section 4.4.4, we compare in a simulation study the numerical behavior of the described methods to estimate EDR directions. We also evaluate the different procedures of selection of the useful components of x. In Section 4.4.5, we apply the most efficient one on real data from a genetic framework. Finally, some concluding remarks are given in Section 4.4.6. 4.4.2 SIR in determined and underdetermined cases 4.4.2.1 Brief review of usual SIR Let β be a p × K matrix defined by β = (β 1 , . . . , β K ). The EDR subspace is thus spanned by β.

Inverse regression step. The basic principle of the SIR method is to reverse the role of y and x, that is, instead of regressing the univariate variable y on the multivariate • We run SIR-QZ for {H 1 , . . . , H N H } = {5, . . . , 15}. In Figure 4.14(c), we plot y against the corresponding estimated indices. This graphic exhibits a structure which is similar to the one in Figure 4.13.

• We finally apply SIR-MP with H = 10. We observe in Figure 4. 14(d), which shows how Y and the indices produced by the SIR-MP method are related, that this method also fails to recover the shape of Figure 4.13.

Thus, Figure 4.14 shows that for this data set, SIR-QZ provides better estimations of the indices than RSIR, SR-SIR and SIR-MP.

To quantify such conclusions, we can use a criterion that measures how X ′ β and X ′ β are close form each other, for a given estimate β. Let X = (x 1 , . . . , xn ) and let P be the The p-dimensional random vector x is generated from the law N (0, Σ) where Σ is a p × p matrix created as follows. Let D = diag(d 1 , . . . , d p ) and for j = 1, . . . , p define d j by

Let also U 1 be a p × p matrix satisfying vec(U 1 ) ∼ N (0, I pp ). See Remark 10 for more information about the vec operator. Let U 2 be the p × p matrix corresponding to an orthogonal basis made of eigenvectors of U 1 U ′ 1 . We then choose Σ = U 2 DU ′ 2 . Model (4.32) and the model from the simulation study in [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF] are the same except that p = 20 in the latter one. For this value of p and n = 500, RSIR outperforms the other methods presented in [START_REF] Zhong | RSIR: regularized sliced inverse regression for motif discovery[END_REF].

Let SIR-I be the method that produces the estimate βSIR-I . We define it as the matrix that contains the K eigenvectors of M which are related with the first K eigenvalues. This is equivalent to a classical SIR estimation, using Σ = I p . We evaluate the performances of RSIR, SR-SIR, SIR-QZ, SIR-MP and SIR-I when trying to accurately estimate X ′ β from samples (X, Y ) of size n = 50 generated from model (4.32). To do so we compute R(m) for each method m and 100 samples. The results are presented in Figure 4.17. There is no approach m such that R(m) appears significantly higher than R(SIR-I). The choice of the parameter s in RSIR and SIR-QZ is thus not sensitive in this case and then comparing these approaches is difficult. That is why, in model (4.29) we introduced a covariance matrix Σ which differs a lot from I p .

Real data application 4.4.5.1 Description of the Dataset

We illustrate our developed approach on a genetic dataset which contains transcriptomic data and genomic data. In this study, we aim at finding genetic causes of variation in the sparsest model, that is with (p 0 , ρ) = (10, 0.9). Figure 4.19 exhibits the selected 11 SNPs for this choice of p 0 and ρ. The threshold (horizontal red line in the figure) is defined as follows: N 1

is the quantile of order (1 -α/2 p ) of the standard normal distribution. It corresponds to the upper bound of the prediction interval of the occurrence of a SNP in the selected model under the hypothesis that none of the SNPs are associated with the gene expression. The level of this interval is fixed at 1α = 0.95 and is corrected by a Bonferroni approach.

On Figure 4.18(b), we plot the dependent variable Hopx versus the index based on the 10 SNPs selected according to our previous comments on Table 4.3. The linear correlation between this index and the one estimated on all the SNPs which is equal to 0.956, highlights the good behaviour of our CCS strategy to select the relevant SNPs. Thus, not surprisingly, we observe the same relation between Hopx gene expression and the estimated indices.

Comparison methods

We compare our approach with three popular multivariate methods for analyzing highdimensionnal datasets: A Lasso approach, the sparse Partial Least Squares (sPLS) and a Bayesian variable selection regression (ESS++). The Lasso method [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF]) often performs poorly in prediction and interpretation especially when n is small and p is large. This technique tends to shrink the regression coefficients towards zero in order to select a sparse subset of covariates and provide a better prediction performance. sPLS [START_REF] Lê Cao | Sparse canonical methods for biological data integration: application to a cross-platform study[END_REF]) seeks for the best linear combination of SNPs to predict the outcome.

To ensure sparsity, sPLS includes a penalty function on some loading coefficients which is equivalent to a restriction on the number of loading vectors and on the number of SNPs, in each vector, that have a non-null coefficient. Both Lasso and sPLS approaches require a preliminary calibration of the tuning parameters which directly affects the number of selected variables, the estimate of the model parameters and therefore the statistical performances of the models. Calibration procedures usually involve the minimization of the mean square error of prediction through V-fold cross validation. In this illustration, we used the leave-one-out crossvalidation method to choose the tuning parameter for both methods. We finally compare our results with ESS++ a Bayesian variable selection approach for linear regression that can analyze single and multiple responses (Bottolo Let consider a sample X = (X 1 , . . . X n ) made of independent and identically distributed random variables generated from the density f . To estimate f from this sample of size n, [START_REF] Parzen | On estimation of probability density function and mode[END_REF] and [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF] introduced the kernel density estimator fK,h , defined for every real t by

where K is called the kernel and is most of time a density function while h is a positive real parameter that controls the smoothness of fK,h . We also refer the interested reader to [START_REF] Scott | Multivariate Density Estimation : Theory, Practice and Visualization[END_REF] and Silverman (1986).

Let N ( fK,h ) be the number of modes of fK,h . To decide how smooth fK,h should be, an approach is to set N ( fK,h ). To do so, one can use the critical bandwidth h crit introduced by [START_REF] Silverman | Using kernel density estimates to investigate multimodality[END_REF] for h. It can be defined by

2 , the bandwidth h crit,k has interesting properties. It can easily be computed and, provided that k ≥ N (f ), allows fK,h to exhibit, in probability, a pointwise convergence, an L 1 -convergence and a uniform convergence toward f . Details can be found in [START_REF] Futschik | On the consistency of kernel density estimates under modality constraints[END_REF], [START_REF] Mammen | Some asymptotics for multimodality tests based on kernel density estimates[END_REF], [START_REF] Devroye | The strong uniform consistency of kernel density estimates[END_REF], [START_REF] Devroye | A Course in Density Estimation[END_REF] and Section 3.2.

These properties are proven using the key point that h crit,k converges in probability toward 0 when K is the Gaussian kernel. To extend them to other kernels, the asymptotic behavior of h crit,k should thus be studied. In this Note, we focus on the case when K is the uniform kernel defined as ∀t ∈ R,

, where 1 is the indicator function. We will prove in the following section that for this kernel, h crit,k does not converge toward 0, after giving some properties about fK,h and N ( fK,h ).

A.1.2 Properties for the uniform kernel

Let us first introduce some additional notations. Let

and

In order to deduce the value of N ( fK,h ), we only need to investigate how the points in A h ∪ B h are ordered because of Proposition 13 below. We write w = card(A h ∪ B h ) ≤ 2n.

We set c h,(i) i∈{1,...,w} as the ordered points in A h ∪ B h . Let us also write b h,(0) = c h,(0) = -∞ and a h,(card(A h )+1) = c h,(w+1) = +∞.

Proposition 13. Let (X 1 , . . . X n ) be a vector of independent random variables generated from f . Let fK,h be the kernel estimator of f for the uniform kernel K. Then, ∀h > 0, ∀i ∈ {0, 1, . . . , w}, the function fK,h is constant on ]c h,(i) , c h,(i+1) [.

The proof is given in Appendix A.3.7.

Remark 16. The reasoning in the proof of Proposition 13 can also be used to obtain the following results:

• ∀i ∈ {1, . . . , w} , c h,(i

Proposition 13 and Remark 16 are illustrated in Figure A.1 where the sample (3, 6.5, 6, 1.5, 5) of size n = 5 is used to compute the function fK,2 . We have here w = 9. Note that N ( fK,2 ) = 3 and that a mode of fK,2 is actually a single point. This mode vanishes if we use a bandwidth slightly less than 2 so that there is a jump in the estimator h → N ( fK,h ) located at h = 2. Figure 16 gives a lead to find N ( fK,h ) from A h and B h , as explained in the following property. Proposition 14. Let (X 1 , . . . X n ) be a vector of independent random variables generated from f . Let fK,h be the kernel estimator of f for the uniform kernel K. The number of modes N ( fK,h ) of fK,h is such that N ( fK,h ) = card (i, j) : a h,(i) ∈]b h,(j-1) , b h,(j) ] and b h,(j) ∈ [a h,(i) , a h,(i+1) [ , where i ∈ {1, . . . , card(A h )} and j ∈ {1, . . . , card(B h )}.

The proof is given in Appendix A.3.8. This characterization of N ( fK,h ) based on the sets A h and B h , together with an argument about totally positive matrices, allows us to show the following theorem:

A.2 Details about SIR-QZ A.2.1 Generalized real Schur decomposition

We work here with real matrices M and Σ. Similarly to the generalized Schur decomposition introduced in Section 4.4.2.3, that produces complex matrices Q and Z, the generalized real Schur decomposition (see Theorem 7.7-2 of [START_REF] Golub | Matrix computations[END_REF]) ensures that we can find Q and Z such that Q M Z is an upper quasi-triangular real matrix and Q ΣZ is an upper triangular real one. An upper quasi-triangular matrix can be defined as the sum of an upper triangular matrix and of a block diagonal matrix where the sizes of the block are 1 × 1 or 2 × 2. For a 1 × 1 diagonal block of the matrix Q M Z, its unique element is called tj if it is located at the jth row and at the jth column of Q M Z. We write tj 1 ,j 2 the element of a 2 × 2 diagonal block, in the j 1 th row and in the j 2 th column of Q M Z. An example of such an upper quasi-triangular matrix is given below:

where * denotes some real values. Let J be made of the elements j ∈ {1, . . . , p} such that tj exists and J c be the set made of j ∈ {1, . . . , p -1} such that tj,j and tj+1,j+1 exist.

For each j ∈ J, let ũj be the element of Q ΣZ at the same location than tj in Q M Z and define similarly ũj 1 ,j 2 for each tj 1 ,j 2 . Thus, we have

det tj,jλũ j,j tj,j+1λũ j,j+1 tj+1,j tj+1,j+1λũ j+1,j+1 .

Hence, for j ∈ J, if ũj = 0, then λ j = tj /ũ j is a real generalized eigenvalue. In addition, for j ∈ J c , [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF] succeeded in finding ( tj , tj+1 ) ∈ C 2 , and (ũ j , ũj+1 ) ∈ R\{0} such that λ j = tj /ũ j and λ j+1 = tj+1 /ũ j+1 are generalized eigenvalues of M and Σ. This leads to vectors t = ( t1 , . . . , tp ) ′ and ũ = (ũ 1 , . . . , ũp ) ′ that are sent back by the QZ algorithm in order to provide generalized eigenvalues.

A.2.2 Algorithm

As explained in Section 4.4.2.3, the QZ algorithm has to be controlled when dealing with singular pencils. The following pseudocode in Scilab language allows the user to do so in the context of sliced inverse regression, for underdetermined cases. -16 , 10 and 10 -10 . The QZ algorithm is implemented in Matlab through the eig function. With Scilab, one should use spec, which is based on the LAPACK library. The R software is able to call functions from this library. Thus the QZ algorithm can be easily tested with this software. Notice that Algorithm A.1 is designed to handle real values of ũj and tj but, as mentioned in Appendix A.2.1, they can be complex. In that case, knowing if the blocks made of tj,j , tj+1,j , tj,j+1 and tj+1,j+1 and of ũj,j , ũj+1,j , ũj,j+1 and ũj+1,j+1 produce unstable eigenvalues is more difficult. As explained in Section 5 of [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF], the QZ algorithm aims at finding stable λ j and λ j+1 corresponding to these 2 × 2 blocks. Because we do not control this procedure, we simply report if the QZ algorithm send back complex values in t. We never encounter this case in the simulation study of Section 4.4.4.

A.2.3 The sliced indices issue

Hereafter, we describe why the Algorithm A.1 produces clustered indices as in Figure

, where every element of the n × n matrix 1 n,n is equal to 1. Notice that X = X Īn and then Σ = 1 n X Īn X ′ . Let S be a n × H matrix made of elements S i,h defined as

We can also write M = X Īn S W S ′ Īn X ′ , where W is defined in Section 4.4.2.2. Because of the structure of S, for any H × α matrix A, SA has at most H distinct rows, so has Īn SA. Let w be the first generalized eigenvector of Īn S W S ′ Īn and Īn associated with the eigenvalue λ. This means that Īn S W S ′ Īn w = λ Īn w which implies that w = Īn w has at most H distinct values and then that w has also at most H distinct values.

Assume that X has full column rank, which is likely to happen when p > n. Then, X +′ w is a generalized eigenvector of M and Σ and X ′ X +′ w = w has at most H distinct values. The eigenvalue that is related to X +′ w is equal to nλ.

If it exists β1 = X +′ w, a generalized eigenvector of M and Σ such that the generalized eigenvalue which is related to β1 is greater than the one corresponding to X +′ w, then we βQZ with H = 10. The horizontal scale was standardized.

But, because w is the first generalized eigenvector of Īn S W S ′ Īn and Īn , it maximizes

over the vectors u of length n, which contradicts the latter equation. Hence, such a β1 does not exist, and the first generalized eigenvector of M and Σ is X +′ w.

In this paragraph, we show that it exists H -1 orthogonal vectors w 1 , . . . , w H-1 such that, for k = 1, . . . , H -1,

which means that the first K generalized eigenvectors of M and Σ are the vectors X +′ w k for k = 1, . . . , K ≤ H -1. We assume n > H. Let sort Y increasingly and reorder the columns of X such that each column corresponds to the appropriate element of Y . This transformation implies that S W S ′ is block diagonal with H blocks. For h = 1, . . . , H, the size of the block h is equal to np h and each element it contains is equal to 1 n 2 ph . Hence, the rank of each block is equal to 1 and each block provides a positive eigenvalue for S W S ′ , which is equal to np h n 2 ph = 1 n . The corresponding eigenvector is made of elements I[y i ∈s h ] √ np h for i = 1, . . . , n. We have now H orthonormal eigenvectors of S W S ′ for the eigenvalue

A.3 Proofs

A.3.1 Proof of Lemma 2

The Markov chain (Z k , S k+1 ) is ergodic and admits a unique invariant measure given by

where U E denotes the uniform distribution on the set E. One may refer the reader to Azaïs et al. for more details on the properties of (Z k , S k+1 ) in the general case. Let 0 ≤ s ≤ t. As a consequence, we have

according to the almost sure ergodic theorem. In particular, the limit for s = t is strictly positive. Besides, since L m (x, •) is a decreasing function, we have R m (x, s) ≤ R m (x, t) at least for m large enough. In this context, we have

A.3.2 Proof of Theorem 9

Recall that, since the bases D, B (1) , . . . , B (q) are assumed to be Σ-orthonormal, we have

Hence, it is well known that the matrix V which maximizes Trace(D ′ ΣBB ′ ΣD) over the set of matrices D such that D ′ ΣD = I K is made of the K generalized eigenvectors j ∈ {j 1 , j 2 } and P B(j) , Σ = P B 2 ,Σ + O P (n -1/2 ), for j ∈ {j 3 , j 4 }. It follows that

Similarly, r( B(j 3 ) , Σ, B(j 4 ) , Σ) tends to 1 in probability when n → ∞.

A .3.6 Proof that (4.28) provides a basis of the EDR space

We have, for k = 1, . . . , K, Σ 1/2 M + Σ 1/2 η k = α k η k . By multiplying by the matrix (M + Σ 1/2 ), we obtain:

Using the definition of b k given in (4.28), we have:

From the comments on the functional SIR context, provided in Section 4.4.2.4, the proof is complete.

A.3.7 Proof of Proposition 13

[ and X (i) i∈{1,...n} be the ordered sequence of the elements of X. We will show that fK,h (u) is neither greater nor lesser than fK,h (v) with a proof by contradiction. Note that for the uniform kernel we have fK,h (u) = h (v), this implies that there exists at least one k ∈ {1, . . . , n}, for which we have X

[ and it is impossible.

A.3.8 Proof of Proposition 14

We will show the equivalence between the presence of a mode between a h,(i) and b h,(j) and the inequality b h,(j-1) < a h,(i) ≤ b h,(j) < a h,(i+1) .

At first, we notice that ordered like this, there is no element of A h or B h that can be between a h,(i) and b h,(j) . This is why the last inequality is equivalent to ∃k ∈ {1, . . . , w -1}, a h,(i) = c h,(k) and b h,(j) = c h,(k+1) , provided that a h,(i) = b h,(j) .

From Proposition 13, fK,h is constant on ]a h,(i) , b h,(j) [, and thanks to Remark 16, it is equivalent to: fK,h is constant on [a h,(i) , b h,(j) ] = [c h,(k) , c h,(k+1) ]. In order for this interval to be a mode, we must prove that there exists ε > 0 for which fK,h is increasing on [c h,(k)ε, c h,(k) [ and decreasing on ]c h,(k+1) , c h,(k+1) + ε], which is also made in Remark 16.

When a h,(i) = b h,(j) = c h,(k) , fK,h is increasing on [c h,(k)ε, c h,(k) [ too and decreasing on ]c h,(k) , c h,(k) + ε]. The mode is reduced to a single point.

where I γ is the γ × γ identity matrix. It is straightforward to show that J is a totally positive matrix since all its minors are not negative.

If ω / ∈ Ω 1 , then D ε,n = D ε,n-1 , because A h and B h stay the same if we build them with (X 1 , . . . , X n ) or with (X 1 , . . . , X n-1 ). Then J = I w and is totally positive.

To conclude, we write ÑK,h : n → ÑK,h (n) = N ( fK,k ). Recall that ÑK,h

2

. Because n → v(D ε,n ) is increasing, ÑK,h is also an increasing function. Let h crit,k,n be the critical bandwidth defined in (A.1) for a sample of size n, then we have: ∀h < h crit,k,n , ÑK,h (n) > k.

Because ÑK,h increases with n, it follows that, ∀η ∈ N, ∀h < h crit,k,n , ÑK,h (n + η) > k.

Thus, with the definition of h crit,k , ∀η ∈ N, h crit,k,n+η ≥ h crit,k,n , and the proof is complete.

List of works