domain: math info sdv

text:
Merci !

Mes premiers remerciements vont à mes trois rapporteurs : Christophe Andrieu, Didier Concordet et Pascal Massart. Je vous suis sincèrement reconnaissante d'avoir accepté de consacrer de votre temps précieux à la lecture de mon manuscrit. Je remercie également tous les autres membres du jury : Florence Forbes, Jean-Michel Marin, Christian Robert, Stéphane Robin et Sophie Schbath. C'est un grand honneur pour moi que vous soyez tous là aujourd'hui. Vous avez toute mon estime et toute ma considération. Merci à vous ! Je remercie ensuite chaleureusement les personnes qui ont contribué de près ou de loin à ce que je prenne la décision d'écrire mon manuscrit en vue d'obtenir le diplôme d'habilitation à diriger les recherches. Merci aussi à celles qui m'ont encouragée et soutenue lors de la rédaction proprement dite et de toutes les étapes qui ont suivi jusqu'à ces derniers jours. Je leur en suis extrêmement reconnaissante. Leurs manisfestations, qu'il s'agisse d'un petit sourire en passant devant mon bureau, d'un petit mot, d'un message ou d'un échange plus profond, furent infiniment précieuses pour moi. J'adresse également mes remerciements à tous mes collaborateurs avec lesquels j'ai pris beaucoup de plaisir à travailler au cours de ces années d'activités de recherche. Sans eux, ce manuscrit ne serait pas tel qu'il est. Merci aussi à tous les collègues avec lesquels j'ai eu l'occasion de partager des connaissances plus ou moins profondes, d'appréhender des interrogations plus ou moins sérieuses, d'échanger des opinions plus ou moins marquées, de confronter des points de vue plus ou moins déterminés. Ces moments de discussion, et l'émulation qui en découle, contribuent considérablement à mes yeux à l'activité de recherche.

Enfin, mes derniers remerciements vont aux membres de tous les laboratoires dans lesquels j'ai exercé avec plaisir et enthousiasme mes activités de recherche et d'enseignement depuis quinze ans : mes débuts en thèse au laboratoire de mathématiques d'Orsay, puis le CEREMADE de l'Université Paris Dauphine et le laboratoire Statistique et Génome de l'Université Evry Val d'Essonne comme ATER, le laboratoire d'analyse, de géométrie et de leurs applications de l'Université Paris Nord en tant que maître de conférence et actuellement l'unité MaIAGE du centre INRA de Jouy-En-Josas comme chargée de recherche. Dans tous ces laboratoires, j'ai trouvé un environnement harmonieux et apaisé, propice au travail. Je souhaite vivement pouvoir poursuivre mon activité dans des conditions aussi agréables.

Introduction en français

Parcours scientifique

Mon stage de DEA effectué sous la direction de Marc Lavielle fut ma première expérience de recherche en statistique. A cette occasion, nous avons collaboré étroitement avec un géophysicien pour résoudre un problème d'estimation en tomographie dans un contexte à données manquantes en implémentant un algorithme numérique. Ce projet initial rassemble à lui seul les composantes principales de mes activités de recherche depuis lors.

J'ai ensuite effectué une thèse en statistique à l'Université Paris Sud sous la direction de Marc Lavielle. J'ai proposé, étudié et mis en oeuvre des algorithmes pour l'estimation par maximum de vraisemblance dans des modèles à variables latentes non linéaires. J'ai effectué des applications en pharmacologie et en traitement du signal en collaboration avec des scientifiques de ces deux domaines. Après mon doctorat, j'ai exercé en tant que maître de conférence à l'Université Paris Nord au Laboratoire d'Analyse, de Géométrie et de leurs applications pendant quatre ans. J'y ai poursuivi mes activités de recherche en statistique et ai développé de nouvelles collaborations motivées par des applications en analyse d'image et en nutrition. Depuis 2009, je suis chargée de recherche à l'INRA au département Mathématiques et Informatique Appliquées (MIA). J'exerce ma fonction au sein de l'équipe DYNENVIE de l'unité MIA de Jouy-en-Josas. Je m'intéresse à de nouvelles applications telles que l'étude de dynamique de population en épidémiologie ou en agronomie, ou encore l'étude de la croissance des plantes en collaboration avec des scientifques de l'INRA et d'autres organismes. J'oriente et développe ma recherche théorique en statistique pour apporter des réponses pratiques aux questions soulevées par ces applications.

Contexte scientifique

Les problématiques statistiques auxquelles je m'intéresse découlent majoritairement de l'analyse de phénomènes complexes dans lesquels interviennent différentes quantités : certaines sont observables, plus ou moins directement, d'autres pas. Un exemple classique est celui de la déconvolution de signal où l'on observe le signal d'intérêt bruité. L'objectif est d'obtenir à partir des seules quantités observées des informations sur les quantités d'intérêt non observées ou partiellement observées. Les outils probabilistes habituellement utilisés dans ce cas sont les modèles à variables latentes. Les quantités d'intérêt non accessibles sont modélisées par des variables aléatoires non observées aussi appelées latentes. Les quantités accessibles sont modélisées par des variables aléatoires pour lesquelles on observe une réalisation. La dépendance entre les observations et les variables latentes est modélisée par la distribution jointe. Un des enjeux principaux pour le statisticien est de caractériser cette distribution à partir des seules observations. Dans la suite, on supposera que cette distribution jointe est paramétrique. L'objectif est alors de fournir un estimateur des paramètres du modèle, ainsi que de la variance de cet estimateur (par exemple pour construire des intervalles de confiance). L'estimateur du maximum de vraisemblance (EMV) est un de ceux considérés classiquement dans ce contexte. Il est obtenu comme solution du problème d'optimisation en le paramètre de la vraisemblance observée, c'est-à-dire la vraisemblance marginale définie comme l'intégrale de la vraisemblance jointe sur les variables latentes. Dans des modèles complexes, cette intégrale n'admet généralement pas de forme analytique et le problème d'optimisation ne peut être résolu par un calcul direct. Dans ce cas, un algorithme d'optimisation peut permettre d'obtenir une approximation numérique de l'EMV.

Un des algorithmes les plus répandus est l'algorithme Expectation Maximization (EM) proposé par [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Il s'agit d'un algorithme déterministe itératif, chaque itération comportant deux étapes. La première consiste à évaluer l'espérance de la log-vraisemblance complète conditionnellement aux observations et à la valeur courante du paramètre (étape E) ; dans la seconde, la valeur du paramètre est mise à jour en maximisant cette quantité (étape M). Pour des modèles exponentiels, cet algorithme converge vers un maximum local de la vraisemblance observée sous des hypothèses générales de régularité du modèle. Cependant, dans de nombreux modèles, l'espérance conditionnelle de la log-vraisemblance complète n'admet pas d'expression analytique et l'algorithme EM ne peut être implémenté. Des algorithmes alternatifs ont été proposés, soit en approximant la vraisemblance, soit en faisant intervenir une étape de simulation des variables latentes. La plupart de ces algorithmes ne possèdent pas de propriétés théoriques de convergence, ou requièrent des hypothèses de convergence peu réalistes, ou encore nécessitent de très longs temps de calcul.

L'algorithme Stochastic Approximation Expectation Maximization (SAEM) proposé par [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] est une version stochastique de l'algorithme EM qui possède d'intéressantes propriétés combinant celles de l'approximation stochastique et celles de l'algorithme EM. L'étape E de l'algorithme EM est remplacée par deux étapes. Dans la première, une réalisation des variables latentes est simulée selon la loi conditionnelle, dans la seconde, cette réalisation est utilisée pour calculer une quantité auxiliaire approximant l'espérance conditionnelle de la log-vraisemblance complète par un schéma d'approximation stochastique. La convergence presque sûre de cet algorithme vers un maximum local de la vraisemblance observée a été établie sous des hypothèses générales de régularité du modèle (cf. [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF]). L'algorithme SAEM est facile à implémenter et nécessite de faibles temps de calcul. Il peut être appliqué à des modèles complexes sous réserve de savoir simuler des réalisations des variables latentes selon la loi conditionnelle aux observations. Cette condition restreint drastiquement le champ des applications possibles. De façon plus générale, l'étape de simulation des variables latentes requiert une attention spécifique. En effet, des comportements numériques atypiques peuvent apparaître lorsque les variables latentes sont de grande dimension, car leur distribution est souvent multimodale dans ce contexte.

Mes thèmes de recherche principaux concernent le développement, l'analyse et l'implémentation de nouveaux algorithmes stochastiques dérivant de l'algorithme déterministe EM et permettant d'obtenir une approximation numérique d'un estimateur des paramètres pour des modèles probabilistes complexes à variables latentes. Les applications font partie intégrante de mes travaux de recherche, les nouvelles applications posant de nouvelles problématiques, à la fois de modélisation et computationnelles, et peuvent elles mêmes conduire à de nouveaux développements mathématiques et algorithmiques. Cette symbiose entre recherche statistique et applications est un élément clé de la dynamique de mes activités de recherche et passe par des collaborations étroites avec des scientifiques d'autres domaines.

Organisation du manuscrit

L'essentiel de mes travaux de recherche est présenté dans ce manuscrit. J'ai choisi de les regrouper en trois parties comme suit.

Dans la première partie, après une brève description du problème de l'estimation dans des modèles à variables latentes et des algorithmes d'optimisation existants, je présente des développements autour de l'algorithme SAEM et leurs analyses. Au cours de ma thèse dirigée par Marc Lavielle, nous avons proposé d'introduire dans l'étape de simulation de l'algorithme SAEM une méthode de Monte Carlo par chaines de Markov (MCMC) permettant ainsi d'appliquer cet algorithme sans avoir besoin de simuler des réalisations des variables latentes selon la loi conditionnelle [A1]. Cet algorithme, noté SAEM-MCMC dans la suite, conserve toutes les propriétés intéressantes de l'algorithme SAEM initial tout en étant facile à mettre en oeuvre. Nous avons démontré sa convergence presque sûre vers un point critique de la vraisemblance observée sous des hypothèses générales de régularité du modèle et sous une hypothèse forte pour la loi des variables latentes, à savoir être à support compact. Par la suite, en collaboration avec Alain Trouvé et Stéphanie Allassonnière, nous avons relâché cette hypothèse en introduisant une étape supplémentaire de troncature à chaque itération de l'algorithme [A3]. Nous avons obtenu le même résultat de convergence pour cet algorithme avec troncature sous des hypothèses générales de régularité seulement. Par ailleurs, l'algorithme proposé est facile à implémenter et rapide. Nous obtenons également un estimateur de la matrice d'information de Fisher observée [A1].

Dans le cas de variables latentes de grande dimension, la performance des méthodes MCMC classiques décroit rapidement lorsque cette dimension augmente. Motivée par des applications en analyse d'images, je me suis intéressée aux techniques de simulation de variables aléatoires de grande dimension. En collaboration avec Stéphanie Allassonnière, nous avons proposé une version anisotrope de l'algorithme Metropolis Adjusted Langevin (AMALA). Nous avons démontré son ergodicité uniforme. Nous l'avons utilisé comme échantillonneur dans l'algorithme SAEM-MCMC et avons prouvé la convergence presque sûre de cet algorithme vers un point critique de la vraisemblance observée, ainsi qu'un théorème de la limite centrale [A12]. En collaboration avec Gersende Fort, Benjamin Jourdain, Tony Lelièvre et Gabriel Stolz, nous avons étudié l'algorithme de Wang Landau, bien adapté pour simuler des variables en grande dimension. Nous avons démontré son ergodicité uniforme ainsi qu'un théorème de la limite centrale [A11]. Nous avons également analysé son comportement via une étude de simulations [A8].

Dans la deuxième partie, je présente mes travaux relatifs à la modélisation et à l'estimation paramétrique dans des modèles à effets mixtes, en particulier dans les modèles déformables en analyse d'image, ainsi que mes travaux liés à des problématiques de tests en régression gaussienne. En collaboration avec Marc Lavielle, nous avons proposé d'utiliser l'algorithme SAEM-MCMC pour obtenir une approximation numérique de l'EMV dans des modèles à effets mixtes. Ces modèles sont fréquemment utilisés pour l'analyse de mesures longitudinales répétées. Lorsque le modèle est linéaire en les effets aléatoires et que l'erreur résiduelle est gaussienne, la vraisemblance est explicite et l'estimateur du maximum de vraisemblance peut être calculé par une optimisation directe. En revanche, hors de ce contexte spécifique, la vraisemblance n'admet pas forcément de forme analytique, rendant le calcul de l'EMV difficile. Il faut alors recourir à des procédures numériques plus complexes, souvent longues en temps de calcul, et pour lesquelles il n'existe pas toujours de propriétés théoriques garantissant leur convergence. L'algorithme SAEM-MCMC est une solution efficace, rapide et convergente, pour calculer l'EMV dans des modèles à effets mixtes, en particulier non linéaires. Nous l'avons appliqué pour des modèles de courbes de croissance et de pharmacocinétique [A2]. En collaboration avec Alain Trouvé et Stéphanie Allassonnière, nous avons appliqué l'algorithme SAEM-MCMC avec l'étape de troncature à un modèle déformable d'images dans un cadre bayésien, motivé par la petite taille des échantillons disponibles en imagerie médicale [A3]. Un modèle déformable permet de représenter un échantillon d'images d'un objet d'intérêt par une image de référence et des déformations géométriques, de sorte que chaque image de l'échantillon est obtenue comme résultat d'une déformation géométrique de l'image de référence, à une petite erreur près. Les déformations géométriques sont les variables latentes du modèle. Nous avons considéré dans un premier temps des petites déformations linéaires et avons utilisé un échantillonneur de type Gibbs hybride comme méthode MCMC pour des variables multivariées de dimension raisonnable. Cette application est à la source des développements liés aux variables latentes de grande dimension présentés dans la première partie, l'algorithme de Gibbs hybride devenant très gourmand en temps de calcul lorsque la dimension des variables latentes augmente. En collaboration avec Stéphanie Allassonnière, nous avons considéré une extension de ce modèle à un modèle de mélanges. Dans de nombreuses applications, il est nécessaire de contraindre le type de déformations géométriques considérées. Par exemple, les déformations difféomorphes n'autorisent pas les recouvrements dans l'image et empêchent des changements de topologies. Un modèle de mélanges permet de modéliser l'échantillon d'images comme étant issu de plusieurs populations ayant des images de référence topologiquement différentes (par exemple avec ou sans "boucle" pour le chiffre manuscrit 2). Dans le cas de mélange de modèles déformables, la mise en oeuvre de l'algorithme SAEM-MCMC donnent lieu à des comportements numériques instables liés aux états absorbants. Nous avons proposé un algorithme stochastique d'estimation spécifique que nous avons mis en oeuvre et pour lequel nous avons établi la convergence vers un point critique de la vraisemblance [A5]. Par la suite, en collaboration avec Stéphanie Allassonnière et Stanley Durrlemann, nous avons appliqué l'algorithme SAEM-MCMC avec l'échantillonneur AMALA au modèle déformable à grandes déformations difféomorphes [A13]. Le nouvel algorithme révèle tout son potentiel dans ce contexte de très grande dimension. Nous avons également développé une extension du modèle permettant d'optimiser la position des points de contrôle de la déformation, simultanément aux autres paramètres du modèle. Nous avons aussi proposé un critère empirique pour sélectionner un nombre optimal de points de contrôle, dans le but d'optimiser la dimension du modèle. Les applications, réalisées en image sur la base US POSTAL et sur des données d'imagerie médicale 2D et 3D, sont regroupées à la fin de cette partie.

En collaboration avec France Mentré et Adeline Samson, nous avons comparé la pertinence de traitements effectués sur deux groupes de patients. Dans ce cadre, j'ai appliqué l'algorithme SAEM-MCMC à des modèles de pharmacologie. Nous avons appliqué l'algorithme SAEM-MCMC pour effectuer l'estimation de paramètres sous l'hypothèse nulle et sous l'alternative et avons ensuite appliqué un test du rapport de vraisemblance. En collaboration avec Mohamed Sahmoudi, Karim Abed-Meraim, Philippe Ciblat et Marc Lavielle, nous avons appliqué l'algorithme SAEM-MCMC pour estimer les paramètres d'un modèle de séparation de sources pour des signaux à queue lourde. Ces travaux détaillés au chapitre 2 de ma thèse de doctorat [T] et dans le proceeding [P2] ne sont pas présentés dans ce manuscrit.

Motivée au départ par la question du choix de modèle pour les modèles à effets mixtes, j'ai proposé en collaboration avec Sylvie Huet un test d'adéquation de modèles pour tester une hypothèse de type linéaire sur l'espérance d'un vecteur gaussien à erreurs bloc corrélées à structure de covariance connue aux paramètres près [A10]. Nous avons démontré que notre procédure de test est asymptotiquement consistante et puissante sur une large classe d'alternatives. Nous avons également proposé une version bootstrap de notre test et montré sa consistance. Nous avons évalué via une étude de simulation les propriétés non asymptotiques de nos procédures sur des échantillons de taille finie et les avons appliquées à un jeu de données de couverture forestière de Galicie. J'ai également effectué un travail d'appui statistique pour la mise en oeuvre d'une procédure de test multiple en régression gaussienne, réalisé en collaboration avec Renaud Rincent et Alain Charcosset, généticiens à l'INRA, pour étudier la puissance de détection de Quantitative Trait Loci le long du génome dans un modèle à effets mixtes linéaire faisant intervenir la matrice d'apparentement. Nous avons étudié les performances de la procédure en fonction de l'estimateur considéré pour cette matrice [A9]. Ce travail n'est pas présenté dans ce manuscrit.

Dans la troisième partie de ce document, j'ai regroupé mes travaux relatifs à la modélisation et l'estimation en analyse de survie. En collaboration avec Charles El-Nouty, j'ai étudié les modèles de fragilité proposés par [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF] qui sont une extension du modèle de Cox et permettent de prendre en compte l'hétérogénéité présente dans les données de survie en introduisant des variables latentes. Notre intérêt pour ces modèles a été motivé par de nombreux échanges avec des chercheurs de l'Unité de Recherche en Epidémiologie Nutritionnelle (UREN) de l'Université Paris 13 qui ont fréquemment besoin d'analyser ce type de jeux de données. Très riches du point de vue de la modélisation, ces modèles sont relativement complexes du point de vue mathématique et l'estimation des paramètres y est souvent difficile. Dans de nombreux cas, les algorithmes existants n'apportent pas de solutions satisfaisantes. Nous avons proposé d'appliquer aux modèles de fragilité l'algorithme SAEM-MCMC pour obtenir une approximation numérique de l'estimateur du maximum de vraisemblance [A6]. Nous avons démontré que sous des hypothèses générales de régularité vérifiées par les modèles de fragilité classiques, l'algorithme SAEM-MCMC était presque sûrement convergent vers un point critique de la vraisemblance observée. Nous avons comparé les performances de cet algorithme à celles d'autres algorithmes sur des données simulées et réelles. Les résultats numériques montrent un net avantage à l'utilisation de l'algorithme SAEM-MCMC, tant du point de vue de la vitesse d'exécution que de celui de la précision.

En collaboration avec Luc Duchateau et Klaartje Goethals, chercheurs en biométrie à la faculté vétérinaire de l'Université de Ghent, nous avons analysé un jeu de données d'épidémie de mastitis en mettant en oeuvre des modèles de fragilité avec des variables de fragilité multivariées possédant différentes structures de corrélation [R1]. Pouvoir considérer ces modèles complexes permet d'étudier le risque de propagation de la maladie en fonction de la position locale de l'infection. L'estimation des paramètres qui était jusqu'alors impossible en un temps raisonnable dans de tels modèles est réalisée par l'algorithme SAEM-MCMC. Nous comparons quatre modèles emboités en utilisant des tests du rapport de vraisemblance. Nous justifions a posteriori la validité à distance finie de ce test par une étude de simulation.

Dans le cadre d'une collaboration avec Catherine Picon-Cochard, biologiste à l'INRA, j'ai effectué un travail d'appui statistique sur l'utilisation du modèle de Cox pour l'analyse de données morphométriques de racines en fonction de covariables environnementales [A7]. Ce travail n'est pas présenté dans ce manuscrit.

Pour finir, mes conclusions et perspectives de recherche sont rassemblées dans la dernière partie du manuscrit.

Introduction

Scientific path

My first experience in statistical research took place during my master's internship supervised by Marc Lavielle (INRIA Saclay Ile de France). We addressed, in close collaboration with a geophysicist an estimation problem in tomography in a missing data setting by implementing a numerical algorithm. This project covered the principal components of my future research activities.

Then, I did a Phd in statistics at University Paris Sud supervised by Marc Lavielle. I have proposed, studied and implemented algorithms for maximum likelihood estimation in non linear latent variables models. I have applied them in pharmacology and signal processing in collaboration with scientists of these fields. After my Phd, I worked as assistant professor at University Paris Nord at the "Laboratoire d'Analyse, de Géométrie et de leurs applications" during four years. I continued my research activities in statistics and developed new collaborations motivated by applications in image processing and nutrition. Since 2009, I am researcher at INRA in the departement "Mathématiques et Informatique Appliquées" (MIA) in the team DYNENVIE of the unit MIA in Jouy-en-Josas. I am interested in new applications as the study of population dynamics in epidemiology or agronomy, and as the study of plant growth in collaboration with scientists of INRA and of others instituts. I position and develop my theoretical research activities with the objective of providing practical answers to the questions raised by these applications.

Scientific context

I am interested in statistical problems raised by the analysis of concrete complex phenomena in which several quantities are involved: some are observed, more or less directly; others are not. A classical example is the deconvolution problem where the signal of interest is observed up to a noise term. The objective is to extract information on the non-observed or partially observed quantity of interest from the observations. The usual probabilistic tools used in such cases are latent variable models. The quantities of interest that are inaccessible are modeled with unobserved random variables, also referred to as latent or missing variables. The observations are modeled with random variables that are observed. The dependence between observations and latent variables is specified by the joint distribution. One of the main goals for the statistician is therefore to characterize this distribution on the basis of observations alone. In this manuscript, I will assume that this joint distribution is parametric. The objective is then to provide an estimate for the model parameter, as well as for the variance of this estimate (for example, to build confidence intervals). The maximum likelihood estimate (MLE) is one of the estimates usually considered in this context. It is obtained by solving the optimization in the parameter of the observed likelihood, i.e., the marginal likelihood defined as the integral of the joint likelihood over the latent variables. In complex models, this integral admits in general no analytical expression and the optimization cannot be directly solved. In such cases, a numerical value of the MLE may be obtained by using an optimization algorithm.

One of the most common algorithms is the Expectation Maximization (EM) algorithm proposed by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. It is a deterministic iterative algorithm, where each iteration is composed of two steps. The first consists in evaluating the conditional expectation of the complete log-likelihood given the current parameter estimate value (step E); in the second one, the parameter value is updated by maximizing this quantity (step M). For models belonging to the curved exponential family, this algorithm converges toward a local maximum of the observed likelihood under general regularity assumptions on the model. However, in numerous models, the conditional expectation of the complete log-likelihood has no analytical expression and the EM algorithm cannot be implemented. Alternative algorithms have been proposed, either based on approximations of the likelihood or by introducing a simulation step of the latent variable. Most of them have no established theoretical convergence property, require unrealistic assumptions for theoretical convergence, or require very long computation times.

The Stochastic Approximation Expectation Maximization (SAEM) algorithm proposed by [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] is a stochastic version of the EM algorithm that has interesting properties, combining those of stochastic approximation procedures and those of the EM algorithm. Step E of the EM algorithm is divided into two steps. In the first one, a realization of the latent variable is drawn from the conditional distribution; in the second one, this realization serves to calculate an auxiliary quantity that approximates the conditional expectation of the complete log-likelihood through a stochastic approximation procedure. The almost sure convergence of this algorithm toward a local maximum of the observed likelihood is established under general regularity assumptions on the model (cf. [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF]). It is easy to implement and requires only short computation times. It can be implemented in complex models as soon as it is possible to draw a realization of the latent variable from the conditional distribution. This condition drastically limits the application possibilities. More generally, the simulation step of the latent variable requires close attention since atypical numerical behavior can result from its implementation, in particular, for multimodal latent variable distributions or high dimensional ones.

My main research subjects deal with the development, the analysis and the implementation of stochastic algorithms derived from the deterministic EM algorithm that makes it possible to calculate a numerical value of a parameter estimate for complex latent variable models. Applications are an integral part of my research activities. New problematics in modeling and in computation arise from applications and raise new theoretical and algorithmic developments. This symbiosis between statistical research and applications is a key element of the dynamics of my research activities and is based on close collaboration with scientists in other fields.

Manuscript organization

In this manuscript, I describe my main contributions, broken down into three parts as follows.

In the first part, after a short description of the estimation problem in general latent variable models and the existing algorithms used to address it, I present developments and analysis of the SAEM algorithm. During my PhD which was supervised by Marc Lavielle, we proposed the introduction of a Monte Carlo Markov Chain (MCMC) method into the simulation step of the SAEM algorithm. This makes it possible to implement the algorithm without being able to draw a realization of the latent variable from the conditional distribution [A1]. This algorithm, designated as SAEM-MCMC below, maintains the good theoretical and practical convergence properties of the SAEM algorithm and, at the same time, extends its application possibilities. We first proved the almost sure convergence of the generated estimate sequence toward the maximum likelihood estimate under usual regularity assumptions on the model and a strong assumption on the latent variable support, i.e., having a compact support. Later, in collaboration with Alain Trouvé (ENS Cachan, CMLA) and Stéphanie Allassionnière (Ecole Polytechnique, CMAP), we relaxed this strong assumption by adding a truncation step to each iteration of the algorithm [A3]. We obtained the same convergence result under general regularity assumptions alone. Moreover, the resulting algorithm is easy to implement and very fast. We also obtained an estimate of the observed Fisher information matrix.

However, considering latent variables of high dimension, the efficiency of classical MCMC methods quickly decreases as the dimension increases. Motivated by complex application settings in image analysis, I focused on sampling methods adapted to this context. In collaboration with Stéphanie Allassionnière, we proposed a new Anisotropic Metropolis Adjusted Langevin Algorithm (AMALA). We proved its uniform ergodicity. We also used it as a sampler in the SAEM-MCMC and proved the almost sure convergence of the generated estimate sequence toward the maximum likelihood estimate under usual regularity assumptions, as well as a central limit theorem [A12]. In addition, in collaboration with Gersende Fort (CNRS, TELECOM ParisTechTelecom, LTCI), Benjamin Jourdain, Tony Lelièvre and Gabriel Stolz, (Ecole des Ponts ParisTech, CERMICS), we studied another sampling algorithm that was well adapted to the high dimension setting, known as the Wang Landau algorithm. We proved its uniform ergodicity as well as a central limit theorem [A11]. We also studied its behavior by performing numerical simulations [A8].

The second part covers my contributions concerning with modeling and parametric estimations in mixed effects models, in particular, in the deformable template model in image analysis, and with the testing problematic in Gaussian regression models.

In collaboration with Marc Lavielle, we proposed using the SAEM-MCMC algorithm to obtain a numerical value of the MLE in mixed effects models. These models are particularly used to analyze repeated longitudinal data. When the model is linear in the random effects and the residual random error Gaussian, the likelihood has an explicit analytical expression and the maximum likelihood estimate can be calculated through some direct optimization. On the other hand, outside this specified context, the likelihood usually does not admit an explicit analytical form, making it difficult to evaluate the MLE. It is then necessary to appeal to more complex numerical procedures, often time-consuming, and not always provided with theoretical convergence properties. The SAEM-MCMC algorithm is an efficient solution, fast and convergent, to evaluate the MLE in mixed effects models, in particular, in non-linear ones. We used it for parameter estimation in growth curve models and in pharmacodynamic models [A2]. In collaboration with Alain Trouvé and Stéphanie Allassonnière,we implemented the SAEM-MCMC algorithm with the additional truncation step to estimate parameters in a deformable template model in a Bayesian setting useful in medical imaging where samples are usually small [A3]. Deformable template models allow us to represent a sample of images with a reference image called template and geometrical deformations, so that each image of the sample is obtained as the result of the geometrical deformation of the template, up to a small error term. The geometrical deformations are the latent variables of such a model. We first considered the case of small linear deformations and used a hybrid Gibbs sampler as the MCMC method. This application motivated the development of high dimensional latent variables presented in the first part of this manuscript, since the hybrid Gibbs sampler becomes very time-consuming as the dimension of the latent variable increases. Later, in collaboration with Stéphanie Allassonnière,we proposed a specific algorithm for parameter estimation in the multicomponent model which is a mixture of deformable template models motivated by a crucial modeling issue [A5]. In numerous applications, it is necessary to constrain the type of geometrical deformations considered, in particular, so that the diffeomorphic deformations do not allow overlapping and prevent topological changes from occuring between the template and the observation, thus creating the need of a mixture model. Since the SAEM-MCMC algorithm is sensitive to numerical phenomena such as trapping states in such a case, we proposed a specific stochastic estimation algorithm and implemented it. We also established its convergence property toward a local maximum of the observed likelihood. Finally, in collaboration with Stéphanie Allassonnière and Stanley Durrlemann, we applied the SAEM-MCMC algorithm provided with the AMALA sampler to the very complex deformable template model using large diffeomorphic deformations [A13]. The new algorithm takes on its full meaning in this high dimensional setting. We also developed an extension that allowed us to optimize the position of the control points of the deformation, simultaneously with the estimation of the other model parameters. We proposed an empirical criterion to select an optimal number of control points as well, leading to the optimization of the model dimension. All the applications related to deformable template models and the different estimation algorithms were performed on the US POSTAL handwritten digit database and 2D and 3D medical image databases. The corresponding experimental results are presented in Section 5.5. I also implemented the SAEM-MCMC algorithm to carry out the parameter estimation in pharmacological models and signal processing. In collaboration with France Mentré and Adeline Samson,we compared the efficiency of treatments given to two patient groups. We used the SAEM-MCMC algorithm to carry out the parameter estimation under the null hypothesis and under the alternative, and subsequently applied the likelihood ratio test. In collaboration with Mohamed Sahmoudi, Karim Abed-Meraim, Philippe Ciblat and Marc Lavielle,we applied the SAEM-MCMC algorithm to carry out the parameter estimation in blind source separation for heavy tailed signals. These contributions are detailed in Chapter 2 of my PhD manuscript [T] and in the proceedings [P2], respectively, and are not presented in this manuscript.

Motivated at the beginning by the model choice issue for mixed effects models, we finally proposed, in collaboration with Sylvie Huet a goodness-of-fit test for testing a linear hypothesis on the expectation of a Gaussian vector with block correlated errors with a known covariance structure up to some parameters [A10]. We established that our test procedure was asymptotically of the nominal level and consistent over a large class of alternatives. We also proposed a bootstrap version of our procedure. Using a simulation study, we evaluated the finite sample size properties of our procedures and applied them to a forest cover dataset of Galicia.

I also did a statistical support study on a multiple test procedure in Gaussian regression, in collaboration with Renaud Rincent and Alain Charcosset, geneticians at INRA, to study the power of detection of Quantitative Trait Loci along the genome in a non linear mixed effects model involving the kinship matrix. We have studied its performances according to the estimate used for the kinship matrix [A9]. This work is not presented in this manuscript.

In the third part, I included my contributions dealing with modeling and estimation in survival analysis. In collaboration with Charles El-Nouty, we considered the frailty models introduced by [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF] which are an extension of the Cox model that takes the heterogeneity that exists in survival data into account by introducing latent variables. We were interested in these models since we have been in contact with researchers at the Unité de Recherche en Epidémiologie Nutritionnelle (UREN) at the University of Paris 13, who often need to analyze this type of dataset. Very rich from a modeling point of view, these models are very complex from a mathematical point of view, and the estimation task is often very difficult. In numerous cases, the existing algorithms do not provide satisfactory solutions. We proposed the application of the SAEM-MCMC algorithm to frailty models to evaluate the MLE [A6]. We proved that under general regularity assumptions fulfilled by classical frailty models, the SAEM-MCMC algorithm is almost surely convergent toward a local maximum of the observed likelihood. We compared the performances of this algorithm with others that exist in the literature on simulated data and on a real set of bladder cancer data. The numerical results showed a net advantage when using the SAEM-MCMC algorithm, both in terms of the accuracy of the limit as well as the computation time.

In collaboration with Luc Duchateau and Klaartje Goethals (Faculty of Veterinary Medicine, Ghent University), we have analyzed a mastitis epidemic dataset using frailty models with a frailty vector of size four with different covariance structures [R1]. Assessing the correlation structure in cow udder quarter infection times allows us to analyze the propagation risk of the disease as a function of the position of the infection. The parameter estimation that could not be performed in a reasonable time in such models until now, is done using the SAEM-MCMC algorithm. We compared four nested models using likelihood ratio tests. Using simulation studies, we justified a posteriori the use of likelihood ratio tests for finite sample size.

I also did a statistical support study in collaboration with Catherine Picon-Cochard on the application of the Cox model for analyzing morphometric plant root data as a function of environmental covariables [A7]. This work is not presented in this manuscript.

Finally, my conclusions and perspectives can be found in the last part of this manuscript.

Part I

Extending the Stochastic Approximation Expectation Maximization Algorithm

The first part is composed of my research dealing with the development, the study and the implementation of a new stochastic version of the Expectation Maximization (EM) algorithm. Related studies include [A1,A3,A8,A11,P1,P2,T,A12].

In the first part, after a short description of the estimation problem in general latent variable models and the existing algorithms used to address it, I present developments and analysis of the SAEM algorithm. During my PhD which was supervised by Marc Lavielle, we proposed the introduction of a Monte Carlo Markov Chain (MCMC) method into the simulation step of the SAEM algorithm. This makes it possible to implement the algorithm without being able to draw a realization of the latent variable from the conditional distribution [A1]. This algorithm, designated as SAEM-MCMC below, maintains the good theoretical and practical convergence properties of the SAEM algorithm and, at the same time, extends its application possibilities. We first proved the almost sure convergence of the generated estimate sequence toward the maximum likelihood estimate under usual regularity assumptions on the model and a strong assumption on the latent variable support, i.e., having a compact support. Later, in collaboration with Alain Trouvé (ENS Cachan, CMLA) and Stéphanie Allassionnière (Ecole Polytechnique, CMAP), we relaxed this strong assumption by adding a truncation step to each iteration of the algorithm [A3]. We obtained the same convergence result under general regularity assumptions alone. Moreover, the resulting algorithm is easy to implement and very fast. We also obtained an estimate of the observed Fisher information matrix.

However, considering latent variables of high dimension, the efficiency of classical MCMC methods quickly decreases as the dimension increases. Motivated by complex application settings in image analysis, I focused on sampling methods adapted to this context. In collaboration with Stéphanie Allassionnière, we proposed a new Anisotropic Metropolis Adjusted Langevin Algorithm (AMALA). We proved its uniform ergodicity. We also used it as a sampler in the SAEM-MCMC and proved the almost sure convergence of the generated estimate sequence toward the maximum likelihood estimate under usual regularity assumptions, as well as a central limit theorem [A12]. In addition, in collaboration with Gersende Fort (CNRS, TELECOM ParisTechTelecom, LTCI), Benjamin Jourdain, Tony Lelièvre and Gabriel Stolz, (Ecole des Ponts ParisTech, CERMICS), we studied another sampling algorithm that was well adapted to the high dimension setting, known as the Wang Landau algorithm. We proved its uniform ergodicity as well as a central limit theorem [A11]. We also studied its behavior by performing numerical simulations [A8].

Estimation in the Latent Variable Models

Latent variable models

We consider a latent variable model defined as follows: the observed variable is denoted by y and the latent variable (also referred to as the missing or hidden variable) by z. The observed variables are related to the latent variables. This link is specified through the joint distribution of the complete variable defined by (y, z). We assume that this joint distribution has a density against a given σ-finite Borelian measure µ which belongs to a parametric family denoted by {f (y, z; θ), θ ∈ Θ}, where the parameter θ takes its value in a subset of R p denoted by Θ. This joint distribution defines the parametric probabilistic latent variable model. The most popular ones are hidden Markov models, hierarchical models and mixed effects models. Such models are widely used in many application fields such as tomography, signal processing, pharmacology, genetics, economics and image analysis.

Given the probabilistic latent variable model, the statistician is interested in estimating the parameter θ from the observed data y, the latent variable z being unobserved. In a frequentist estimation approach, a very popular estimator with good asymptotical properties in many classical models is the Maximum Likelihood Estimator (MLE) of the parameter θ, namely the value θ of θ that maximizes the observed likelihood g defined by: g(y; θ) = f (y, z; θ)µ(dz).

(1.1)

We can also consider a Bayesian estimation approach, for example if only few observations are available or if a priori information is available and has to be introduced into the probabilistic model. The parameter θ is considered as a random variable that follows a given prior distribution denoted by q. A useful estimator would then be the Maximum A Posteriori (MAP) estimator, namely the value θ of θ that maximizes the posterior distribution denoted by p and defined by: p(θ|y) ∝ g(y; θ)q(θ).

In the following, I consider a frequentist estimation approach and focus on the Maximum Likelihood Estimator (MLE). However all the considerations developed in this section can be directly transposed to the Bayesian context by replacing the MLE with the MAP estimator, as will be done in Section 5.

As soon as the joint distribution of the complete variable (y, z) is complex, it is often not possible to directly compute the numerical value of such estimators as the solution of an optimization problem. It is then necessary to use numerical algorithms to evaluate an estimated value of this estimator.

Numerical algorithms for estimation

A useful tool to compute a numerical value of the Maximum Likelihood Estimator θ is the Expectation Maximization (EM) algorithm introduced by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. This is an iterative algorithm that generates a sequence of estimated values (θ k ) k converging toward a stationary point of the observed likelihood g under some regularity assumptions of the model (see [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], [START_REF] Vaida | Parameter convergence for EM and MM algorithms[END_REF], [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF]). The heuristic of this algorithm is based on the following strategy: the variable z being unobserved, instead of maximizing the observed log-likelihood log g(y; θ) in θ, it is instead possible to maximize the expectation of the complete log-likelihood conditional to the observed variable y equal to E(log f (y, z; θ)|y; θ), which might be a realistic approximation of log g(y; θ) denoted by l(θ) for a given observation y in the sequel. The two steps of iteration k of the algorithm consist in alternatively updating, first, the conditional expectation:

Q(θ|θ k-1 ) = E(log f (y, z; θ)|y; θ k-1 ) (1.2)
using the current parameter value θ k-1 and, second, the parameter by maximizing Q according to θ k = argmax θ Q(θ|θ k-1 ). The initial value θ 0 is chosen arbitrarily. Applying Jensen's inequality shows that the sequence (log g(y, θ k )) k is non-decreasing where the sequence (θ k ) k is generated by the algorithm. The structural hypothesis always assumed to study its theoretical convergence is that the complete likelihood belongs to the curved exponential family, meaning that the complete likelihood f (y, z; θ) can be written as:

f (y, z; θ) = exp [-ψ(θ) + S(z), φ(θ) ] ,
where •, • denotes the Euclidean scalar product, the sufficient statistics S is a function on R l , taking its values in a subset S of R m and ψ, φ are two functions on Θ (note that S may also depend on y, but we omit this dependency for the sake of simplicity). This condition is usual within the framework of EM algorithm applications and is fulfilled by a large range of complex models. Moreover, [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF] exhibits some strong assumptions that ensure the convergence of the sequence (θ k ) k toward a local maximum of the observed log-likelihood log g.

Nevertheless, one or both steps of this algorithm may not always be easily carried out in complex latent variable models. Thus, alternative algorithms were developed to allow wider range of applications for the EM algorithm.

When the M-step cannot be done analytically, alternative satisfactory solutions have been proposed, e.g., based on the Newton algorithm (see [START_REF] Lange | A gradient algorithm locally equivalent to the EM algorithm[END_REF]) or on conditional maximization (see [START_REF] Meng | On the rate of convergence of the ECM algorithm[END_REF]).

On the contrary, overcoming the difficulty of an untractable E-step leads to alternative solutions that are not as satisfactory. When the computation of the conditional expectation defined in Equation (1.2) is not feasible, two types of alternatives, in particular, have been proposed. The first ones are based on an approximation of the underlying likelihood function, e.g., by its first or second order development. To the best of my knowledge, there is no theoretical convergence property established in that case, probably since it is difficult to control the propagation of the error made by such an approximation. Several authors proposed using a Laplace approximation for the computation of the conditional expectation defined in Equation (1.2). In the case of non-linear mixed effects models, it was proposed by [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF] but very restrictive assumptions were required to prove the theoretical convergence of the generated sequence (θ k ). The second ones are based on the simulation of the unobserved variable. Let us mention the stochastic EM algorithm proposed by [START_REF] Celeux | L'algorithme SEM: Un algorithme d'apprentissage probabiliste pour la reconnaissance de mélange de densités. (The SEM algorithm: An algorithm of probabilistic learning for the determination of mixtures of densities)[END_REF]. The authors proved the almost sure convergence in mean of the generated sequence toward a stationary point of the observed likelihood. They extended their approach to mixture models by adding a simulated annealing step and proved the almost sure convergence toward a local maximum likelihood. Another algorithm was proposed by [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the Poor's Man's data augmentation algorithms[END_REF] who used a Monte Carlo sum to approximate the conditional expectation defined in Equation (1.2) at each iteration leading to the so-called Monte Carlo EM (MCEM) algorithm. [START_REF] Fort | Convergence of the Monte Carlo Expectation Maximization for curved exponential families[END_REF] established its almost sure convergence toward a local maximum of the observed likelihood but this algorithm is time consuming since it is necessary to simulate a huge quantity of auxiliary variables at each iteration.

The Stochastic Approximation Expectation Maximization (SAEM) algorithm

The Stochastic Approximation Expectation Maximization (SAEM) algorithm was proposed by [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] and combined the idea of simulating the unobserved variable to the nice convergence property of the stochastic approximation. It is based on the construction of a sequence (Q k (θ)) k , which will asymptotically produce a good approximation of the conditional expectation defined in Equation (1.2). The E-step of the EM algorithm is divided into two steps. Consider the iteration k. First, a realization z k of the latent variable is sampled from the conditional distribution π θ k-1 (z|y) of z conditional to y using the current value of the parameter estimate θ k-1 . Second the quantity Q k is updated through a stochastic approximation procedure using the realization z k . Let (γ k ) k be a decreasing positive step size sequence. The algorithm remains to:

Initialization step: Initialize θ 0 in a fixed compact set.

Then, for all k ≥ 1 the k th iteration consists in three steps :

Simulation step: simulate z k from the conditional distribution π θ k-1 (•|y). Stochastic approximation step: compute the quantity

Q k (θ) = Q k-1 (θ) + γ k-1 (log f (y, z k ; θ) -Q k-1 (θ)),
(1.3)

Maximization step: update the parameter value according to θ k = argmax θ Q k (θ). [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] proved the almost sure convergence of the sequence (θ k ) k toward a local maximum of the observed likelihood by using martingale theory results. This algorithm therefore has very interesting convergence properties and requires only short computation times. Nevertheless, it suffers from several drawbacks. Most of them are common to other algorithms of this type. First, since the limit is a local maximum, it may depend on the initialization point, making it necessary to run the algorithm several times with different initializations. Second, the nature of the limit point has to be investigated to be sure that it is a global maximum. Finally, the last drawback, specific to this algorithm, is the requirement of simulating the unobserved variable from the conditional distribution, significantly limiting its application.

Using the Monte Carlo Markov Chain method in the SAEM algorithm

During my PhD, which was supervised by Marc Lavielle, we proposed a new version of the classical SAEM algorithm that does not require the simulation of the unobserved variable from the conditional distribution.

I provide the details of the Monte Carlo Markov Chain method below and its insertion into the SAEM algorithm.

Monte Carlo Markov Chain (MCMC) method

Given some objective distribution π to be reached, a MCMC procedure consists in generating an ergodic Markov chain having it as limiting distribution (see [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF], [START_REF] Robert | Méthodes de Monte Carlo par chaînes de Markov[END_REF]). The most common type of MCMC algorithm is the Metropolis Hastings algorithm. Its transition probability consists in sampling a candidate from a proposal (also referred to as instrumental) distribution q and accepting it with a probability equal to the acceptance ratio defined as:

α(z, z ) = min π(z |y)q(z , z) π(z|y)q(z, z ) , 1 (2.1)
This provides a transition kernel Π of this form: for any Borel set A ∈ B:

Π(x, A) = A α(x, z)q(x, z)dz + 1 A (x)
X

(1 -α(x, z))q(x, z)dz .

(2.2)

This procedure always accepts the new value z when the likelihood ratio π(z |y) q(z,z ) is larger than the previous one.

Another well-known MCMC algorithm for multivariate random variables is the Gibbs sampler (see [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF], [START_REF] Robert | Méthodes de Monte Carlo par chaînes de Markov[END_REF]). Its transition probability consists in updating each coordinate at a time by simulating it from its conditional distribution, conditional to all of the other coordinates. This algorithm is particularly well adapted to the Bayesian context when using conjugate laws, since simulation tasks are carried out directly. Otherwise, it is possible to use a Metropolis Hastings step in one simulation step of the Gibbs sampler, leading to the so-called hybrid Gibbs sampler. However, since the dimension of the variable becomes huge, the computation time drastically increases and numerical problems such as trapping states may appear. This often limits its implementation, as will be seen in Section 5.2.

The SAEM-MCMC algorithm

Let us assume that we are not able to simulate the missing data z k at iteration k from the conditional distribution π θ k (•|y). In such a case, we propose using a MCMC procedure to simulate the value z k of the latent variable at iteration k. In practice only the simulation step (S-step) of the SAEM algorithm is modified; the other two remain unchanged. In fact, instead of simulating the unobserved variable from the conditional distribution, we simulate it from the transition kernel Π θ of an ergodic Markov Chain that has this conditional distribution π θ as stationary distribution. Since we consider only parametric models P in this manuscript that belong to the curved exponential family, the stochastic approximation can be made either on the complete log-likelihood or on the sufficient statistics S of the model using a positive step-size sequence (γ k ) k∈N . Finally, we update the parameter in the M-step.

Thus the algorithm can be summarized as follows:

Initialization step: Initialize θ 0 , s 0 and z 0 in fixed compact sets.

Then, for all k ≥ 1 the k th iteration consists in three steps :

Simulation step: simulate z k from the transition probability Π θ k-1 (z k-1 , •).
Stochastic approximation step: compute the quantity

s k = s k-1 + γ k-1 [S(z k ) -s k-1 ],
(2.3)

Maximization step: update the parameter value according to:

θ k = θ(s k ) where θ(s) = argmax θ [-ψ(θ) + S(z), φ(θ) ] . (2.4)
This algorithm still has interesting theoretical properties and its implementation is fast and easy. Moreover, it can be applied to many complex latent variable models. Assuming that the latent variable z has a compact support, we prove the almost sure convergence of the sequence (θ k ) k toward a stationary point of the observed likelihood under some usual regularity assumptions on the model, on the transition kernel of the MCMC method, ensuring its ergodicity, and on the step size sequence (see [A1]).

Theorem 2.1. Assume that some usual regularity assumptions on the model, the transition kernel and the step size sequence are fulfilled. Assume in addition the assumption (C): the sequence (s k ) k≥0 takes its values in a compact subset of S. Then, w.p. 1, lim k→+∞ d(θ k , L) = 0 where d(x, A) denotes the distance of x to the closed subset A and L = {θ ∈ Θ, ∂ θ l(y; θ) = 0} is the set of stationary points of l.

We apply a result of [START_REF] Benveniste | On the consistency of Fréchet means in deformable models for curve and image analysis[END_REF] to control small random Markovian perturbations. To relax the restrictive assumption of bounded support for the latent variable, we proposed, in collaboration with Stéphanie Allassonnière and Alain Trouvé, to add a truncation step to the algorithm (see [A3]). In fact, to prove the convergence of the sequence generated through the stochastic approximation procedure, it is first necessary to prove its stability i.e. it stays in some given compact set. The assumption of bounded support for the latent variable allowed us to prove this in the algorithm without using the truncation procedure. Introducing a more general technique that involves truncation on random boundaries allows us to prove this without assuming that the support of the latent variable distribution is bounded. Our proof is based on the general stability and convergence results for stochastic algorithms with truncation on random boundaries given in [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF]. The main technical point is that in the presence of unbounded latent variables, the usual regularity conditions as a function of the parameters for the solutions of the Poisson equations for the Markovian dynamic cannot be verified and have to be relaxed.

The truncation on random boundaries can be formalized as follows. Let (K q ) q≥0 be a sequence of increasing compact subsets of S so that ∪ q≥0 K q = S and K q ⊂ int(K q+1 ), for all q ≥ 0. Let ε = (ε k ) k≥0 be a monotone non-increasing sequence of positive numbers and K a compact subset of R N . We construct a sequence ((z k , s k )) k≥0 as follows. As long as the stochastic approximation does not wander outside the current compact set and is not too far from its previous value, we run the SAEM-MCMC algorithm. As soon as one of these conditions is not satisfied, we reinitialize the sequences of z and s using a projection, we increase the size of the compact set and continue the iterations until convergence (for more details see [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF]). This is detailed in the following steps:

Initialization step: Initialize z 0 and s 0 in two fixed compact sets.

Then, for all k ≥ 1 the k th iteration consists in four steps :

Simulation step: Draw one new element z of the non-homogeneous Markov Chain with respect to the kernel with the current parameters Π θ k-1 and starting at z k-1 .

z ∼ Π θ k-1 (z k-1 , •) .
Stochastic approximation step: . Compute

s = s k-1 + γ ζ k-1 (S(z) -s k-1
) .

(2.5) [2005].

Maximization step: Update the parameters using (2.4).

The index κ denotes the current active truncation set, the index ζ is the current index in the sequences γ and ε and the index ν denotes the number of iterations since the last projection.

Estimation of the observed Fisher Information Matrix

An estimation procedure should generate a point estimate θ together with the covariance of the estimate (e.g., to enable construction of confidence sets for the true parameter value). The asymptotic theory for maximum likelihood estimation, when established, ensures that:

√ n( θ -θ ) → n→∞ N (0, I(θ ) -1 ) , (2.6)
where θ is the MLE and I(θ ) is the observed Fisher Information Matrix. Thus, an estimate of the asymptotic covariance of θ is the inverse of the observed Fisher information matrix -∂ 2 θ l( θ).

Thanks to the maximum likelihood estimator obtained with the SAEM-MCMC algorithm, it is possible to simultaneously obtain an estimation of the Fisher Information Matrix. In [A1] we propose a method to estimate this matrix by using the fact that the gradient (the Fisher score function) and the Hessian (observed Fisher Information) of the log-likelihood l can be almost directly obtained from the simulated missing data z. Using the so-called Fisher identity, the Jacobian of the log-likelihood of the observed data l(θ) is equal to the conditional expectation of the complete data likelihood:

∂ θ l(θ) E[∂ θ log f (y, z; θ)|y; θ]
where ∂ θ denotes the differential with respect to θ. By analogy with the implementation of the SAEM algorithm, the following approximation scheme is proposed:

∆ k = ∆ k-1 + γ k ∂ θ log f (y, z (k) ; θ k ) -∆ k-1 .
Using Louis' missing information principle [START_REF] Louis | Finding the observed information matrix when using the EM algorithm[END_REF]), the Hessian of l at θ is the observed Fisher Information matrix ∂ 2 θ l(θ) that may be expressed as:

∂ 2 θ l(θ) = E θ [∂ 2 θ log f (y, z; θ)] + Cov θ [∂ θ log f (y, z; θ)].
where Cov

θ (ψ(z)) E θ [(ψ(z) -E θ (ψ(z)))(ψ(z) -E θ (ψ(z))) t ].
Using this expression, it is possible to derive the following stochastic approximation procedure to approximate ∂ 2 θ l(θ):

G k = G k-1 +γ k ∂ 2 θ log f (y, z (k) ; θ k )+∂ θ log f (y, z (k) ; θ k )∂ θ log f (y, z (k) ; θ k ) t -G k-1 H k = G k -∆ k ∆ k t .
Knowing that the algorithm proposed above converges to a limiting value θ and that l is regular enough, (-H k ) converges to the inverse of the observed Fisher Information Matrix -∂ 2 θ l(θ ) (see [A1]).

Toward an efficient sampling step in high dimension

A new Anisotropic Metropolis Adjusted Langevin Algorithm (AMALA)

In collaboration with Stéphanie Allassonnière, we proposed a new anisotropic version of the well-known Metropolis Adjusted Langevin Algorithm (MALA) [A12].

Let us first recall the steps of the Metropolis Adjusted Langevin Algorithm (MALA), which is a particular case of the Metropolis Hastings algorithm (see [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF]). The focus is on optimizing the proposal distribution.

Let us denote by π the pdf of the target distribution with respect to the Lebesgue measure on X , an open subset of R l . We assume that π is positive continuously differentiable. At each iteration k of this algorithm, a candidate X c is simulated with respect to the Gaussian distribution with expectation X k + σ 2 2 D(X k ) and covariance σ 2 Id l where X k is the current value,

D(x) = b max(b, |∇ log π(x)|) ∇ log π(x) , (3.1)
Id l is the identity matrix in R l and b > 0 is a fixed truncation threshold. Note that the truncation of the drift D was already suggested in [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF] to provide more stability.

The Gaussian proposal of the MALA algorithm is optimized with respect to its expectation guided by the Langevin diffusion. One step further is to optimize also its covariance matrix. A first step in this direction was proposed in [START_REF] Atchadé | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF]. The covariance matrix of the proposal is given by a projection of a stochastic approximation of the empirical covariance matrix. It produces an adaptive Markov chain. This process involves some additional tuning parameters that have to be calibrated. Since our goal is to use this sampler in an estimation algorithm, the sampler has a different target distribution (depending on the current estimate of the parameter) at each iteration. Therefore, the optimal tuning parameter may be different along the iterations of the estimation process. Although we agree with the idea of using adaptive chains, we prefer taking advantage of the dynamic of the estimation algorithm. For these reasons, we propose a sampler in the spirit of [START_REF] Atchadé | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF], [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF] or [START_REF] Marshall | An adaptive approach to langevin MCMC[END_REF] without providing an adaptive chain. The adaption will result from the dependency of the target distribution with respect to the parameters of the model that are updated throughout the estimation algorithm. The proposal remains a Gaussian distribution, but both the drift and the covariance matrix depend on the gradient of the target distribution. At the k th iteration, we are provided with X k . The candidate is sampled from the Gaussian distribution with expectation X k + δD(X k ), and the covariance matrix δΣ(X k ) denoted in the sequel N (X k + δD(X k ), δΣ(X k )) where Σ(x) is given by:

Σ(x) = εId l + D(x)D(x) T , (3.2)
D is defined in Equation (3.1) and ε > 0 is a small regularization parameter. Note that the threshold parameter b leads to a symmetric positive definite covariance matrix with bounded non zero eigenvalues.

We introduce the gradient of log π into the covariance matrix to provide an anisotropic covariance matrix depending on the amplitude of the drift at the current value. When the drift is large, the candidate is likely to be far from the current value. This big step may not be of the right amplitude and a large variance will allow for more flexibility. Moreover, this makes it possible to explore a larger area around these candidates, which would not be possible with a fixed variance. On the other hand, when the drift is small in a particular direction, it means that the current value is within a region of high probability for the next value of the Markov chain. Therefore, the candidate should not move too far neither with a large drift nor with a large variance. This makes it possible to extensively sample around large modes, which is of particular interest. This covariance also makes it possible to treat the directions of interest with different amplitudes of variance, as is already the case with the drift. It also provides dependencies between coordinates since the directions of large variances are likely to be different from the Euclidean axis. This is taken into account here by introducing the Gram matrix of the drift into the covariance matrix.

Since our purpose was to plug the AMALA sampler into the SAEM-MCMC algorithm, we have to consider a parametric family of target densities (π s ) s and the corresponding transition kernels (Π s ) s . In particular, we exhibit new assumptions that enable us to prove that the AMALA sampler is uniformly geometrically ergodic when considering a parametric family of target densities (π s ) s and the corresponding transition kernels (Π s ) s .

We require a usual assumption on the stationary distributions namely the so-called super-exponential property given by: (B1) For all s ∈ S, the density π s is positive with continuous first derivative such that:

lim |x|→∞ n(x).∇ log π s (x) = -∞ (3.3) and lim sup |x|→∞ n(x).m s (x) < 0 (3.4)
where ∇ is the gradient operator in R l , n(x) = x |x| is the unit vector pointing in the direction of x and m s (x) = ∇πs(x) |∇πs(x)| is the unit vector in the direction of the gradient of the stationary distribution at point x.

We assume also some regularity properties of the stationary distributions with respect to s.

(B2) For all x ∈ X , the functions s → π s and s → ∇ x log π s are continuous on S.

We now define for some

β ∈]0, 1[, V s (x) = c s π s (x) -β where c s is a constant so that V s (x) ≥ 1 for all x ∈ X . Let also V 1 (x) = inf s∈S V s (x) and V 2 (x) = sup s∈S V s (x).
Let us assume conditions on V 2 : (B3) There exists b 0 > 0 such that, for all s ∈ S and x ∈ X , V b0 2 is integrable against Π s (x, .) and

lim sup b→0 sup s∈S,x∈X Π s V b 2 (x) = 1 . (3.5)
We obtain the following result:

Proposition 3.1. Assume (B1-B3). Let K a compact subset of S. There exist a function V ≥ 1, a set C ⊆ X , a probability measure ν such that ν(C) > 0 and there exist constants

λ ∈]0, 1[, b ∈ [0, ∞[ and ε ∈]0, 1] such that for all s ∈ K : Π s (x, A) ≥ εν(A) ∀x ∈ C ∀A Borel set , (3.6) Π s V (x) ≤ λV (x) + b1 C (x) . (3.7)
The proof is performed in three steps. We first prove the existence of small sets being any compact subset of R l . Then, we prove the Drift condition for each transition kernel Π s with a function V s for all s ∈ S following the lines of [START_REF] Jarner | Geometric ergodicity of Metropolis algorithms[END_REF] and [START_REF] Atchadé | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF]. The fact that both the drift and the covariance matrix are bounded even depending on the gradient of log π s enables partially similar proofs.

Finally, the most technical step consists in exhibiting a single function V built from the family of functions {V s } s satisfaying the Drift condition for all kernels Π s for s ∈ S.

We highlight the efficiency of the AMALA sampler by comparing its mixing properties with those of the MALA sampler. We used both algorithms to sample from a 100 dimensional normal distribution with a zero mean and a non-diagonal covariance matrix. Its eigenvalues range from 1 to 10. The eigendirections are chosen randomly. Ten examples of autocorrelations of both chains are plotted in Figure 1 where we can see that there is a significant benefit of using the anisotropic sampler. To evaluate the weight of the anisotropic term D(x)D(x) T in the covariance matrix, we compute its amplitude (as its non zero eigenvalue since it is a rank one matrix). We see that it is of the same order as the diagonal part on average and increases up to 15 times more. This shows the importance of the anisotropic term.

Coupling the AMALA sampler and the SAEM algorithm

In collaboration with Stéphanie Allassonnière, we proposed using the AMALA sampler in the simulation step of the SAEM-MCMC algorithm presented in Section 2.2 [A12]. Thus, at each iteration k of the algorithm, simulated values of the missing data are drawn from the transition probability of the AMALA algorithm described in Section 3.1 with the current value of the parameters. The others steps remain unchanged.

We proved that the parameter estimate sequence generated by the AMALA-SAEM algorithm converges almost surely toward a stationary point of the likelihood under some regularity assumptions.

Theorem 3.1. Assume some regularity conditions on the model and some usual conditions on the step size sequences. Assume that the family of conditional density probability functions {π θ(s) (•|y), s ∈ S} satisfies (B1-B3).

Let K be a compact subset of X and K 0 a compact of S. Then, for all z 0 ∈ K and s 0 ∈ K 0 , we have lim k→∞ d(θ k , L) = 0 a.s. where (θ k ) k is the sequence generated by the AMALA-SAEM Algorithm and

L {θ ∈ Θ, ∂ θ l(θ) = 0}.
The proof follows the same lines as the one given in [A3]. In particular, we first prove the sufficient usual Drift conditions (cf. [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF]). We also prove in details that the transition kernel Π s is Lipschitz in s ∈ S. This proof extends the one proposed in [START_REF] Andrieu | On the ergodicity properties of some adaptive MCMC algorithms[END_REF] to kernels and stationary distributions both depending on the parameter s ∈ S. The complete algorithm involves only three parameters: b, the threshold for the gradient that appears in the expectation as well as in the covariance matrix; δ, the scale on this gradient; and ε, a small regularization parameter to ensure a positive definite covariance matrix. The scale δ can be easily optimized in terms of the data we are dealing with to adapt to the range of the drift. The value of the threshold b is, in practice, never reached. We also established a Central Limit Theorem for the parameter estimate sequence generated by the AMALA-SAEM algorithm.

Theorem 3.1 ensures that the number of re-initializations of the sequence of stochastic approximation of the AMALA-SAEM Algorithm is finite almost surely. We can therefore consider only the non truncated sequence when we are interested in its asymptotic behavior. Moreover there are a priori multiple possible limiting points so we need to restrict our attention to the set of trajectories that converge to a given limiting point θ * = θ(s * ).

Let us introduce some usual assumptions in the spirit of these of [START_REF] Delyon | Stochastic approximation with decreasing gain: convergence and asymptotic theory[END_REF].

(N1) The function h is C 1 in some neighborhood of s * with first derivatives Lipschitz and J the Jacobean matrix of the mean field h in s * has all its eigenvalues with negative real part. (N2) Let g θ(s) be a solution of the Poisson equation g -Π θ(s) g = H s -p θ(s) (H s ) for any s ∈ S. There exists a bounded function w such that

w -Π θ(s * ) w = g θ(s * ) g T θ(s * ) -Π θ(s * ) g θ(s * ) (Π θ(s * ) g θ(s * ) ) T -U (3.8)
where the deterministic matrix U is given by :

U = E θ(s * ) g θ(s * ) (z)g θ(s * ) (z) T -Π θ(s * ) g θ(s * ) (z)Π θ(s * ) g θ(s * ) (z) T .
(3.9) (N3) The step size sequence (γ k ) is decreasing and satisfies γ k = 1/k α with 2/3 < α < 1.

Theorem 3.2. Under the assumptions of Theorem 3.1 and under (N1)-(N3), the sequence (s k -s * )/ √ γ k converges in distribution to a Gaussian random vector with zero mean and covariance matrix Γ where Γ is the solution of the following Lyapunov equation:

U + JΓ + ΓJ T = 0.
Moreover, denoting θ * = θ(s * ), we have:

1 √ γ k (θ k -θ * ) → L N (0, ∂ s θ(s * )Γ∂ s θ(s * ) T ).
This proof follows the lines of the proof of Theorem 25 of [START_REF] Delyon | Stochastic approximation with decreasing gain: convergence and asymptotic theory[END_REF]. However several of its assumptions are not satisfied by our stochastic approximation. Therefore we exhibit lighter assumptions leading to the same final result when combined with Drift and Hölder conditions.

Convergence study of the Wang Landau Algorithm

In collaboration with Gersende Fort, Benjamin Jourdain, Tony Lelièvre and Gabriel Stoltz, we study the Wang Landau algorithm, from a theoretical point of view as well as from a practical one [A11,A8].

I was interested in this study, having in mind to use the Wang Landau algorithm as sampler in the SAEM-MCMC algorithm. However I do not investigate further this possibility for two main reasons. First, the AMALA sampler developed in Section 3.1 gave very satisfaying theoretical and practical results when used as sampler in the SAEM-MCMC algorithm; second, using adaptive sampler such as the MALA version proposed by [START_REF] Atchadé | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF] into the estimation algorithm let to less performant results. One possible interpretation is that the dynamics of the parameter update in the estimation algorithm plays the role of the adaption process in the sampling step. Therefore I do not consider the coupling of the Wang Landau algorithm with the SAEM algorithm.

Introduction

The Wang-Landau algorithm belongs to the class of free energy biasing techniques (see [START_REF] Lelièvre | Computation of free energy profiles with adaptive parallel dynamics[END_REF]) which have been introduced in computational statistical physics to efficiently sample thermodynamic ensembles and to compute free energy differences. These algorithms can be seen as adaptive importance sampling techniques, the biasing factor being adapted on-the-fly in order to flatten the target probability measure along a given direction. Let us explain this with more details.

Let π be a multimodal probability measure over a high-dimensional space X ⊆ R D . Classical algorithms to sample π such as a Metropolis-Hastings procedure with local proposal moves typically converge very slowly to equilibrium since high probability regions are separated by low probability regions. Averages have to be taken over very long trajectories in order to visit all the modes of the target probability measure π.

The idea of free energy biasing techniques is to flatten the target probability along a well-chosen direction through an importance sampling procedure in order to more easily sample π. More precisely, assume that we are given a measurable function O defined on X and with values in a low dimensional compact space, or in a discrete space. Let us introduce O * π the image of the measure π by O: for any test function ϕ on the image O(X) of X by O, O(X) ϕ(y)O * π(dy) = X ϕ(O(x))π(dx). The free energy biased probability measure π is defined by the two following properties:

(i) the image O * π of π by O is the uniform measure on O(X).

(ii) for each y ∈ O(X), the conditional distributions of x given O(x) = y under π(dx) and π (dx) coincide i.e. there exists a measurable function h

: O(X) → R + such that π (dx) = h(O(x))π(dx).
The bottom line of free energy biasing techniques is that it should be easier to sample π than to sample π since, by construction, O * π is the uniform probability measure. Then, sampling from π could be obtained by importance sampling from π . The fact that π is indeed much easier to sample than π actually depends on the choice of O. It is not an easy task to define and to design in practice a good choice for O and we do not discuss further these aspects here. This is related to the choice of a "good" reaction coordinate in the physics literature, which is a very debatable subject. We refer for example to [START_REF] Chopin | Free energy methods for Bayesian inference: efficient exploration of univariate Gaussian mixture posteriors[END_REF] for such an analysis in the context of free energy biasing techniques used to sample posterior distributions in Bayesian statistics.

Of course, the difficulty is that in general, O * π is unknown so that it is not possible to sample from π . The idea is then to approximate O * π on the fly in order to, in the longtime limit, sample from π . This is the adaptive feature of these algorithms: the importance sampling factor is computed as time goes, in order to penalize states (namely level sets of O) which have already been visited. To approximate π at a given time, one could either use the occupation measure of the Markov chain up to the current time or one could use an approximation over many Markov chains running in parallel (see [START_REF] Lelièvre | Computation of free energy profiles with adaptive parallel dynamics[END_REF], [START_REF] Minoukadeh | Potential of mean force calculations: a multiple-walker adaptive biasing force approach[END_REF]).

In terms of mathematical analysis, approximations based on many replicas in parallel are typically easier to analyze, since they can be related in the limit of infinitely many replicas to mean field models for which powerful longtime convergence analysis techniques can be used. We refer for example to [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method[END_REF] and [START_REF] Lelièvre | Long-time convergence of an adaptive biasing force method: the bichannel case[END_REF] for such an analysis. The convergence analysis and, more importantly, the study of the efficiency of free energy biased techniques for approximations based on the occupation measure are much more involved since correlations in time of the Markov process play a crucial role. Our aim is to propose a convergence analysis for the Wang-Landau algorithm.

Description of the Wang-Landau algorithm

The Wang-Landau algorithm both computes a penalty sequence {θ n , n ≥ 0} approximating in the longtime limit the probability measure O * π and samples draws {X n , n ≥ 0} distributed in the longtime limit according to π . The update of the penalty sequence follows a Stochastic Approximation algorithm (see [START_REF] Benveniste | On the consistency of Fréchet means in deformable models for curve and image analysis[END_REF], [START_REF] Robbins | A stochastic approximation method[END_REF]) and is of the form

θ n+1 = θ n + γ n+1 H n (X n+1 , θ n ) .
Different strategies about the field H n and the adaption schedule {γ n , n ≥ 1} have been proposed in the literature. In the original paper of [START_REF] Wang | Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram[END_REF], the authors came up with a stochastic adaption schedule hereafter called flat histogram Wang-Landau. In this procedure, the updating parameter γ n remains constant up to the random time when the sampling along the chosen order parameter O is approximately uniform, the "amount of uniformity" being measured according to the current value of γ n . Then γ n is lowered and a new updating procedure of the weights starts with a constant stepsize. Another strategy consists in a deterministic update of the adaption sequence {γ n , n ≥ 1}.

We now describe the Wang landau algorithm we studied. Let us consider a partition X 1 , . . . , X d of X in d ≥ 2 elements, and define, for any i ∈ {1, . . . , d},

θ (i) def = Xi π(x)λ(dx) .
(3.10)

In the following, X i will be called the i-th stratum. Each weight θ (i), which is assumed to be positive, gives the relative likelihood of the stratum X i ⊂ X. In practice, the partitioning could be obtained by considering some smooth function ξ : X → [a, b] and defining, for i = 1, . . . , d -1,

X i = ξ -1 [α i-1 , α i ) , (3.11)
and

X d = ξ -1 ([α d-1 , α d ]), with a = α 0 < α 1 < . . . α d = b (possibly, a = -∞ and/or b = +∞).
In our previous notation, the order parameter is thus the discrete function O defined by ∀x ∈ X, O(x) = i if and only if x ∈ X i .

(3.12)

In the following we consider a function O with values in a discrete finite set {1, . . . , d}. Then, we have

π (dx) = 1 d d i=1 1 O(x)=i θ (i) π(dx) , (3.13)
where

θ (i) = π ({x ∈ X, O(x) = i}) = O * π(i) for i ∈ {1, . . . , d}.
The above discussion motivates the fact that the weights θ (i) typically span several orders of magnitude, some sets X i having very large weights, and other ones being very unlikely under π. Besides, trajectories bridging two very likely states may need to go through unlikely regions. To efficiently explore the configuration space, and sample numerous configurations in all the strata X i , it is therefore a natural idea to resort to importance sampling strategies and reweight appropriately each subset X i . A possible way to do so is the following. Let Θ be the subset of (non-degenerate) probability measures on {1, . . . , d} given by

Θ = θ = (θ(1), . . . , θ(d)) 0 < θ(i) < 1 for all i ∈ {1, . . . , d} and d i=1 θ(i) = 1 .
For any θ ∈ Θ, we define the probability density π θ on (X, X ) (endowed with the reference measure λ) as

π θ (x) = d i=1 θ (i) θ(i) -1 d i=1 π(x) θ(i) 1 Xi (x) . (3.14)
This measure is such that the weight of the set X i under π θ is proportional to θ (i)/θ(i). In particular, all the strata X i have the same weight under π θ . Unfortunately, θ is unknown and sampling under π θ is typically unfeasible. The Wang-Landau algorithm precisely is a way to overcome these difficulties: at each iteration of the algorithm, a weight vector θ n = (θ n (1), . . . , θ n (d)) is updated based on the past behavior of the algorithm and a point is drawn from a Markov kernel P θn with invariant density π θn . The intuition for the convergence of this algorithm is that if {θ n , n ≥ 0} converges to θ then the draws are asymptotically distributed according to the density π θ . Conversely, if the draws are under π θ , then the update of {θ n , n ≥ 0} is chosen such that it converges to θ .

We now describe precisely the algorithm we considered. Let {γ n , n ≥ 1} be a [0, 1)-valued deterministic sequence. For any θ ∈ Θ, denote by P θ a Markov transition kernel onto (X, X ) with unique stationary distribution π θ (x)λ(dx); for example, P θ is one step of a Metropolis-Hastings algorithm with target probability measure π θ (x)λ(dx).

Consider an initial value X 0 ∈ X and an initial set of weights θ 0 ∈ Θ (typically, in absence of any prior information, θ 0 (i) = 1/d). Define the process {(X n , θ n ), n ≥ 0} as follows: given the current value (X n , θ n ),

• Draw X n+1 under the conditional distribution P θn (X n , •);

• Set i = O(X n+1 )
where O is given by (3.12). The weights are then updated as

θ n+1 (i) = θ n (i) + γ n+1 θ n (i) (1 -θ n (i)) , θ n+1 (k) = θ n (k) -γ n+1 θ n (k) θ n (i) for k = i. (3.15)
Note that since γ n ∈ [0, 1), θ n ∈ Θ for any n ≥ 0. The update of the probability vector θ n can be recast equivalently into the stochastic approximation framework upon writing

θ n+1 = θ n + γ n+1 H(X n+1 , θ n ) , (3.16)
where

H : X × Θ → [-1, 1] d is defined componentwise by H i (x, θ) = θ(i) (1 Xi (x) -θ(O(x))) .
(3.17)

The updating strategy (3.15) (or equivalently (3.16)) is a modification of the original Wang-Landau algorithm obtained by (i) using a deterministic schedule for the evolution of the step-sizes used to modify the values of the weights (instead of reducing the value of these step-sizes at random times when the empirical frequencies of the strata are sufficiently uniform: this is the flat histogram version of the Wang-Landau algorithm) and (ii) linearizing at first order in γ n the update of the weight θ n .

Concerning this second point, the standard Wang-Landau update is

θ n+1 (i) = θ n (i) 1 + γ n+1 1 Xi (X n+1 ) 1 + γ n+1 θ n (O(X n+1 ))
.

(3.18)

The update (3.15) is obtained from (3.18) in the limit of small γ n .

Convergence of the Wang-Landau algorithm

Despite the Wang-Landau algorithm has been successfully applied for many problems of practical interest, there are many open questions about its longtime behavior and its efficiency. Such a longtime behavior study relies on the convergence of stochastic approximation algorithms with Markovian inputs (see [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF], [START_REF] Benveniste | On the consistency of Fréchet means in deformable models for curve and image analysis[END_REF]) combined with the convergence of adaptive Markov chain Monte Carlo samplers [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF]; for both parts, the stability of the sequence {θ n , n ≥ 0} is a fundamental property. Stability here means that the sequence {θ n , n ≥ 0} remains in a compact subset of the probability measures on {1, . . . , d} with support equal to the support of O * π.

We consider here the Wang-Landau algorithm with a deterministic adaption sequence {γ n , n ≥ 1} for a precise definition of the algorithm) and address both the convergence of {θ n , n ≥ 0} to O * π and the convergence of {X n , n ≥ 0} to π . More precisely, we prove first that the sequence {θ n , n ≥ 0} is stable, which is a crucial point for applications: no ad hoc stabilization techniques (such as truncation at randomly varying bounds [START_REF] Chen | Convergence and robustness of the Robbins-Monro algorithm truncated at randomly varying bounds[END_REF]) is required. We also prove the almost-sure convergence of {θ n , n ≥ 0} as well as a Central Limit Theorem. We then prove the ergodicity and a strong law of large numbers for the draws {X n , n ≥ 0}.

We adopt this linear update (3.15) for the stability and the convergence analysis. The main advantage is that it makes the proof of convergence simpler; nevertheless, since γ n converges to zero, these stability and convergence results are unchanged and could be proved along the same lines for the standard .

The proof of the convergence of the Wang-Landau algorithm relies on its reformulation (3.16) as a stochastic approximation procedure. Since the draws {X n , n ≥ 1} satisfy for any measurable non-negative function f :

E [f (X n+1 )|F n ] = P θn f (X n ) , (3.19)
where F n denotes the σ-field σ(θ 0 , X 0 , X 1 , . . . , X n ).

The main difficulty, when proving the almost-sure convergence of such algorithms, is the stability, namely how to ensure that the sequence {θ n , n ≥ 0} remains in a compact subset of Θ. We use a traditional approach to answer this question: we first prove that our algorithm satisfies a recurrence property i.e. the sequence {θ n , n ≥ 0} visits infinitely often a compact subset of Θ; we then show that there exists a Lyapunov function with respect to the mean-field function h

: Θ → [-1, 1] d h(θ) = X H(x, θ) π θ (x) λ(dx) =   d j=1 θ (j) θ(j)   -1 (θ -θ) , (3.20)
with strong enough properties so that the recurrence property implies stability. Different strategies based on truncations are proposed in the literature to circumvent the stability problem (see e.g. [START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF]). The most popular technique is the truncation to a fixed compact set but this is not a satisfactory solution since the choice of this compact is delicate: a necessary condition for convergence is that the compact contains the unknown desired limit. An adaptive truncation has been proposed by [START_REF] Chen | Convergence and robustness of the Robbins-Monro algorithm truncated at randomly varying bounds[END_REF] which avoids the main drawbacks of the deterministic truncation approach.

We first prove that, under conditions on the target density π and the step-size sequence {γ n , n ≥ 1}, the algorithm (3.16) is recurrent, so that such truncation techniques are not required.

We detailed here the assumptions on the Metropolis dynamics and on the adaption rate. Our conditions fall into three categories: conditions on the equilibrium measure (see A1), on the transition kernels {P θ , θ ∈ Θ} (see A2) and conditions on the step-size sequence {γ n , n ≥ 1} (see A3). It is assumed that A1 The probability density π with respect to the measure λ is such that 0 < inf X π ≤ sup X π < ∞. In addition, inf 1≤i≤d θ (i) > 0 where θ is given by (3.10).

The minorization condition on π certainly is the most restrictive assumption: it is introduced in order to prove the recurrence of the algorithm (3.16). This condition can be removed by adding a stabilization step to (3.16) (such as a truncation technique at random varying bounds [START_REF] Chen | Convergence and robustness of the Robbins-Monro algorithm truncated at randomly varying bounds[END_REF], [START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF]) in order to ensure the recurrence.

A2 For any θ ∈ Θ, P θ is a Metropolis-Hastings transition kernel with invariant distribution π θ dλ, where π θ is given by (3.14), and with symmetric proposal kernel q(x, y)λ(dy) satisfying inf X 2 q > 0.

The minorization condition on q implies that the transition kernels {P θ , θ ∈ Θ} are uniformly (geometrically) ergodic.This property allows a simple presentation of the main ingredients for the limiting behavior analysis of the algorithm. Extensions to a more general case could be done by using the same tools as in [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF] (see also [Andrieu et al., 2005, Section 3]) and controlling the dependence upon θ of the ergodic behavior.

A3 The sequence {γ n , n ≥ 1} is a [0, 1)-valued deterministic sequence such that 1. {γ n , n ≥ 1} is a non-increasing sequence and lim n γ n = 0; 2. n γ n = ∞; 3. n γ 2 n < ∞.
We first proved the following result.

Proposition 3.2. Under A1 and A2, there exists ρ ∈ (0, 1) such that for all θ ∈ Θ, for all x ∈ X and for all A ∈ X , it holds:

P θ (x, A) ≥ ρ A π θ (x) λ(dx) , (3.21) sup θ∈Θ sup x∈X P n θ (x, •) -π θ dλ TV ≤ 2(1 -ρ) n , (3.22)
where for a signed measure µ, the total variation norm is defined as

µ TV = sup {f : sup X |f |≤1} |µ(f )| .
We state that, almost surely, there exists a compact subset of Θ such that θ n belongs to this compact subset for infinitely many n. For any n ≥ 0, set

θ n = min 1≤j≤d θ n (j) .
(3.23)

We prove the following theorem:

Theorem 3.1. Assume A1, A2 and A31. Then,

P lim sup n→∞ θ n > 0 = 1 . (3.24)
The proof is based on the following consideration. The value of the smallest weight increases when the chain goes into the corresponding stratum (see the updating formula (3.15)). Under the stated assumptions, we prove that the chain {X n , n ≥ 0} returns in the strata of smallest weights often enough for the smallest weight to remain isolated from 0.

Then we addressed the almost-sure convergence of the sequence {θ n , n ≥ 0} to θ . Theorem 3.2. Assume A1, A2 and A3. Then, P lim n→∞

θ n = θ = 1.
The proof relies on [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF] which provides sufficient conditions for convergence of stochastic approximation techniques. The first step consists in rewriting the weight update (3.16) as

θ n+1 = θ n + γ n+1 h(θ n ) + γ n+1 H(X n+1 , θ n ) -h(θ n ) , (3.25)
where h is given by (3.20). The heuristic idea is that, if the step-size quickly is sufficiently small, and the Metropolis dynamics converges sufficiently fast to equilibrium for θ fixed (a result given by Proposition 3.2), the update of θ n is indeed close to an update with the averaged drift h(θ n ). However, in order for the updates of the weights to be non-negligible, the step-sizes should not be too small. The balance between these two opposite effects is encoded in the conditions A32-3.

From a technical viewpoint, the proof of the theorem relies on two main tools. The first one is to show that the function V : Θ → R + given by

V (θ) def = d i=1 θ (i) log θ (i) θ(i) (3.26)
is a Lyapunov function with respect to the mean-field h, namely ∇V (θ), h(θ) < 0 for θ = θ and ∇V (θ ), h(θ ) = 0 (here, •, • denotes the scalar product in R d ). This motivates the fact that {θ n , n ≥ 0} may converge to θ . The second important result establishes that the remainder term γ n+1 (H(X n+1 , θ n ) -h(θ n )) in (3.25) vanishes in some sense. This step is quite technical and requires regularity-in-θ of the transition kernels P θ and the invariant distributions π θ . The conclusion then follows from [Andrieu et al., 2005, Theorem 2.3] and Theorem 3.1.

Finally we studied the asymptotic behavior of the chain {X k , k ≥ 0}. We established the following result.

Theorem 3.3. Assume A1, A2 and A3. Then, for any bounded measurable function f ,

lim n→∞ E [f (X n )] = X f (x) π θ (x) λ(dx) , (3.27) 1 n n k=1 f (X k ) a.s. -→ X f (x) π θ (x) λ(dx) . (3.28)
This theorem shows that the distribution of the sample X n converges to π θ (x)λ(dx), where, we recall

π θ (x) = 1 d d i=1 π(x) θ (i) 1 Xi (x) .
Moreover, the empirical mean of the samples {f (X k ), k ≥ 0} converges to f π θ dλ. Hence, although the weights θ n evolve in the adaptive algorithm, ergodic averages can be thought of as averages with fixed weights θ .

In many practical cases, averages with respect to π are of interest. In this case, the Wang-Landau procedure is used as some adaptive importance sampling strategy. In order to obtain averages according to π along a trajectory of the algorithm, some reweighting has to be considered. A natural strategy is to use some stratified-type weighted sum of the samples {X k , k ≥ 1}:

I n (f ) def = d d i=1 θ n (i) 1 n n k=1 f (X k )1 Xi (X k ) .
We also prove the following result:

Theorem 3.4. Assume A1, A2 and A3. Then for any bounded measurable function f ,

lim n→∞ d E d i=1 θ n (i) f (X n ) 1 Xi (X n ) = X f (x) π(x) λ(dx) , (3.29) I n (f ) a.s. -→ X f (x) π(x) λ(dx) . (3.30)
There are of course many other reweighting strategies. We have discussed only one possible choice. May be the above estimator is not the best one.

We also state a Central Limit Theorem on the error (θ n -θ ). We show that the rate of convergence depends upon the step-size sequence {γ n , n ≥ 1} and discuss an averaging strategy in order to reach the optimal rate of convergence. An additional assumption is required on the sequence {γ n , n ≥ 1}:

A4 lim n γ n √ n = 0
, and one of the following condition holds:

1. log(γ n /γ n+1 ) = o(γ n ); 2. log(γ n /γ n+1 ) ∼ γ n /γ with γ > d/2.
Theorem 3.5. Assume that A1, A2, A3 and A4 hold. Then {γ

-1/2 n (θ n -θ ) , n ≥ 1} converges in distribution to a centered Gaussian distribution with variance-covariance matrix σ 2 U where σ 2 = d/2 in case A4(1) and σ 2 = γ d/(2γ -d) in case A4(2), U def = X H θ (x) H T θ (x) -P θ H θ (x) P θ H T θ (x) π θ (x) λ(dx) , (3.31) and H θ def = n≥0 P n θ (I -π θ ) H(•, θ ) = n≥0 P n θ (H(•, θ ) -h(θ )) .
Notice that H θ is the Poisson solution associated to the pair (P θ , H(•, θ )), namely H θ is a solution to: find g : X → R such that

g -P θ g = H(•, θ ) - X H(x, θ ) π θ (x) λ(dx) .
By Proposition 3.2 and the results of [Meyn and Tweedie, 2009, Chapter 17], such a function exists and is unique up to an additive constant.

Theorem 3.5 shows that the rate of convergence depends upon the step-size sequence {γ n , n ≥ 1}: when γ n = γ /n α for α ∈ (1/2, 1], the maximal rate of convergence is reached with α = 1 and the rate is O(n -1/2 ). When γ n = γ /n, one could be interested in optimizing the variance-covariance matrix: introducing a gain matrix Γ in the algorithm (3.16) yields the update θn+1 = θn + γ n+1 Γ H(X n+1 , θn ) .

It is proved in [Benveniste et al., 1990, Proposition 4 p.112] that for a large family of gain matrix (socalled "admissible gains") a Central Limit Theorem still holds for the sequence of random variables { √ n(θ n -θ ), n ≥ 0}, the minimal variance-covariance is equal to d 2 U and is reached with Γ = dγ -1 Id. Since Γ is a scalar matrix, this discussion evidences that the minimal variance-covariance matrix

d 2 U is reached when choosing γ n = d/n.
From a practical point of view, it is known that stochastic approximation algorithms are more efficient when the step-size sequence decreases at a slow rate: in the polynomial schedule, this means that γ n = γ /n α with α close to 1/2. As shown by Theorem 3.5, this yields a slower rate of convergence. Nevertheless, combining Wang-Landau update with an averaging technique allows to reach the optimal rate of convergence and the optimal variance-covariance matrix: by applying [Fort, 2014, Theorem 1.4

], it can be proved that { √ n 1 n n k=1 θ k -θ , n ≥ 1} converges in distribution to a centered Gaussian distribution with variance-covariance matrix d 2 U .

Numerical studies

We complete our theoretical convergence results by a simulation study to discuss the efficiency of the Wang-Landau procedure.

This algorithm is actually known to be useful in metastable situations, namely when the original Markov chain (with transition kernel P θ0 ) remains trapped for very long times in some regions (called the metastable states). Metastability is one of the major bottleneck of standard Markov Chain Monte Carlo techniques, since ergodic averages should be considered over very long times in order to obtain accurate results. Our aim is to show that in such a metastable situation, the Wang-Landau algorithm indeed is an efficient sampling procedure. We consider a toy model composed of only three strata: two large probability strata (the metastable states) separated by a low probability stratum (the transition state). We analyze theoretically the first exit times out of a metastable state.

We show that the Wang-Landau algorithm allows to rapidly escape from a metastable state, namely from a large probability stratum surrounded by small probability strata. We are able to precisely quantify the time the system needs to go from the first metastable state to the second one. We show in particular that the exit time is dramatically reduced with the Wang-Landau dynamics compared to the corresponding non-adaptive dynamics.

Finally we show that (most of) the results obtained for the very simple three-state model are still valid for a less simple example inspired by target measures used in computational statistical physics.

Part II

Inference in mixed effects models, in deformable template models and in

Gaussian regressions

The second part included my contributions dealing with modeling and parametric estimation in mixed effects models, in particular in the deformable template model in image analysis, and with the testing problematic in Gaussian regression models. Related studies include [A2,A3,A4,A5,A10,A12,A13,P3,P4,P5,T].

In collaboration with Marc Lavielle, we proposed using the SAEM-MCMC algorithm to obtain a numerical value of the MLE in mixed effects models. These models are particularly used to analyze repeated longitudinal data. When the model is linear in the random effects and the residual random error Gaussian, the likelihood has an explicit analytical expression and the maximum likelihood estimate can be calculated through some direct optimization. On the other hand, outside this specified context, the likelihood usually does not admit an explicit analytical form, making it difficult to evaluate the MLE. It is then necessary to appeal to more complex numerical procedures, often time-consuming, and not always provided with theoretical convergence properties. The SAEM-MCMC algorithm is an efficient solution, fast and convergent, to evaluate the MLE in mixed effects models, in particular, in non-linear ones. We used it for parameter estimation in growth curve models and in pharmacodynamic models [A2]. In collaboration with Alain Trouvé and Stéphanie Allassonnière,we implemented the SAEM-MCMC algorithm with the additional truncation step to estimate parameters in a deformable template model in a Bayesian setting useful in medical imaging where samples are usually small [A3]. Deformable template models allow us to represent a sample of images with a reference image called template and geometrical deformations, so that each image of the sample is obtained as the result of the geometrical deformation of the template, up to a small error term. The geometrical deformations are the latent variables of such a model. We first considered the case of small linear deformations and used a hybrid Gibbs sampler as the MCMC method. This application motivated the development on high dimensional latent variables presented in the first part of this manuscript, since the hybrid Gibbs sampler becomes very time-consuming as the dimension of the latent variable increases. Later, in collaboration with Stéphanie Allassonnière,we proposed a specific algorithm for parameter estimation in the multicomponent model which is a mixture of deformable template models motivated by a crucial modeling issue [A5]. In numerous applications, it is necessary to constrain the type of geometrical deformations considered, in particular, so that the diffeomorphic deformations do not allow overlapping and prevent topological changes from occuring between the template and the observation, thus creating the need of a mixture model. Since the SAEM-MCMC algorithm is sensitive to numerical phenomena such as trapping states in such a case, we proposed a specific stochastic estimation algorithm and implemented it. We also established its convergence property toward a local maximum of the observed likelihood. Finally, in collaboration with Stéphanie Allassonnière and Stanley Durrlemann, we applied the SAEM-MCMC algorithm provided with the AMALA sampler to the very complex deformable template model using large diffeomorphic deformations [A12,A13]. The new algorithm takes on its full meaning in this high dimensional setting. We also developed an extension that allowed us to optimize the position of the control points of the deformation, simultaneously with the estimation of the other model parameters. We proposed an empirical criterion to select an optimal number of control points as well, leading to the optimization of the model dimension. All the applications related to deformable template models and the different estimation algorithms were performed on the US POSTAL handwritten digit database and 2D and 3D medical image databases. The corresponding experimental results are presented in Section 5.5.

Motivated at the beginning by the model choice issue for mixed effects models, we finally proposed, in collaboration with Sylvie Huet a goodness-of-fit test for testing a linear hypothesis on the expectation of a Gaussian vector with block correlated errors with a known covariance structure up to some parameters [A10]. We established that our test procedure was asymptotically of the nominal level and consistent over a large class of alternatives. We also proposed a bootstrap version of our procedure. Using a simulation study, we evaluated the finite sample size properties of our procedures and applied them to a forest cover dataset of Galicia.

Maximum Likelihood Estimation in Mixed Effects Models

Mixed effects models

Mixed effects models were introduced mainly for modeling repeated longitudinal data (see [START_REF] Davidian | Nonlinear models for Repeated Measurement Data[END_REF]). Such models allow us to analyze responses of a population of individuals that share a global behavior with the exception of some individual variations. Some of them are shared by all the individuals of the population, whereas the others are random, depending on the individuals or possibly on sub-groups of the population. Thus, the model has two types of parameters: global parameters defined as the fixed effects, and parameters that vary among the population defined as random effects. These kinds of observations are usually the result of repeated measurements of some individuals that are repeatedly observed under different experimental conditions. Such settings are very common in practice, for example, in the fields of pharmacokinetics, biological growth, epidemiology and econometry. We refer to [START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF] for more details on mixed effects models.

Let us consider here the following general mixed effects model defined as:

y ij = g(t ij , φ i , β) + h(t ij , φ i , β)ε ij , 1 ≤ i ≤ n , 1 ≤ j ≤ m i (4.1)
where y ij is the jth observation of the ith individual, at some known instant t ij . Here, n is the number of individuals and m i is the number of observations of individual i. The random effects (φ i ) are assumed to be independent identically distributed. We denote by η the parameter of their common distribution. The vector β also denotes unknown population parameters, that do not appear in the distribution of the random effects φ i . The within-group errors (ε ij ) are assumed to be i.i.d. Gaussian random variables with a mean zero and unknown variance σ 2 . We assume that the (ε ij ) and the (φ i ) are mutually independent. The model is said to be nonlinear if the functions g or h are nonlinear in the random effects φ i . Note that mixed effects models are a particular type of latent variable model introduced in Section 1 since the responses (y ij ) are observed and the random effects (φ i ) are unobserved.

Maximum Likelihood Estimation

We choose a frequentist approach and consider the maximum likelihood estimator of the unknown parameter vector θ = (β, η, σ 2 ).

The theoretical study of the MLE has been addressed by [START_REF] Nie | Strong consistency of MLE in nonlinear mixed-effects models with large cluster size[END_REF]. It is a difficult task. The authors studied different asymptotics, as the number of individuals and/or the number of observations per individual tend to infinity. Under some regularity assumptions on the model, they proved the consistency and the asymptotic normality of the MLE.

Let us now focus on the computation of the MLE in practice. For linear mixed effects models, the estimation of the unknown parameters can be treated with the usual EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]), or with a Newton-Raphson algorithm [START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF]). However, nonlinear functions are often more suitable for modeling considerations. Estimating the parameters by maximizing the observed likelihood then requires a specific approach. Different methods, based generally on linearization of the log-likelihood, were suggested for dealing with nonlinear models. A Laplace approximation was proposed by [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF], and a Bayesian approach was proposed by Racine-Poon [1985], [START_REF] Wakefield | Bayesian analysis of linear and nonlinear population models by using the Gibbs sampler[END_REF], [START_REF] Wakefield | The Bayesian analysis of population pharmacokinetic models[END_REF]. [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF] uses a Monte-Carlo EM algorithm, whereas a simulated pseudo maximum likelihood estimator for these specific models was developed by [START_REF] Concordet | A simulated pseudo-maximum likelihood estimator for nonlinear mixed models[END_REF]. These methods are either not proven to be convergent or are very time-consuming. Therefore, applying the SAEM-MCMC algorithm presented in Section 2 is a powerful tool for the estimation task in mixed effects models.

Applications to a growth curve model and to a pharmacodynamic model

The orange tree dataset

We consider the example of orange trees to illustrate our algorithm. This data was studied by [START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF] and is available on S-plus. The data consist in seven measurements of the trunk circumference of each of five orange trees. [START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF] uses a logistic curve to model the trunk circumference y ij of tree i at age x j :

y ij = g(x j , φ i ; β 1 , β 2 ) + ε ij 1 ≤ i ≤ n , 1 ≤ j ≤ m, (4.2) g(x j , φ i ; β 1 , β 2 ) = φ i 1 + exp - xj -β1 β2 . (4.3)
We assume here that the error terms ε ij are independent Gaussian centered variables of variance σ 2 . On the one hand, the asymptotic trunk circumference φ i is treated as a random effect and is assumed to be Gaussian with mean µ and variance τ 2 . On the other hand, the age at which the tree attains half of its asymptotic trunk circumference β 1 and the growth scale β 2 are treated as two fixed effects. Setting

g 1 (φ i ) = φ i and g 2 (β 1 , β 2 , x j ) = 1 1 + exp - xj -β1 β2 , (4.4)
the likelihood of the complete model has the form:

f (y, φ; θ) = (2πσ 2 ) -nm 2 (2πτ 2 ) -m 2 × exp   - 1 2σ 2 i,j (y ij -g 1 (φ i )g 2 (β 1 , β 2 , x j )) 2 - 1 2τ 2 i (φ i -µ) 2   where θ = (β 1 , β 2 , µ, τ 2 , σ 2 ).
Note that the EM algorithm can be implemented in this model. Thus, the value obtained with this algorithm may be considered as the maximum likelihood estimate of θ. The estimation of the parameters after 100 and 1000 iterations with EM and SAEM are displayed in Table 1.

We can observe that EM has almost converged after 100 iterations. In this example, the step size sequence (γ k ) used for SAEM was: γ k = 1 for 1 ≤ k ≤ 100 and γ k = (k -99) -1 for k ≥ 100. After some iterations, the SAEM algorithm has converged to a neighborhood of the MLE of θ. Since γ k = 1 during the first iterations, no further stochastic approximation is performed. Thus, the behavior of the sequence (θ SAEM k ) remains quite perturbed until iteration 100. After that, the introduction of a decreasing step size allows the almost sure convergence of the sequence (θ SAEM k

) to θ M LE .
The Fisher Information of the MLE can also be estimated by using the stochastic approximation scheme presented in Section 2.3. In Table 2, we present the estimated standard deviation of each component of (θ EM ) and (θ SAEM ), obtained after 100 and 1000 iterations. We observe once again that the SAEM algorithm provides a good estimation in just a few iterations.

y ij = φ i 1 + exp - xj -β1 β2
(1 + ε ij ). (4.5)

Note that we are outside the scope of the exponential model. One solution consists in regarding the fixed parameters (β 1 , β 2 ) as the realization of a Gaussian random vector of mean (µ 1 , µ 2 ) and a diagonal covariance matrix with diagonal terms (τ 2 1 , τ 2 2 ). As before, φ = (φ i ) is a sequence of i.i.d. Gaussian random variables of mean µ and variance τ 2 .

It is important to observe that we do not change the model by doing this. The fixed effects remain fixed effects, since we still consider only one vector (β 1 , β 2 ) for the whole population.

Comparisons with other methods on a pharmacodynamic model

In this section, we consider the nonlinear population pharmacodynamic model used by [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF] for comparing the MLEs obtained with the EM algorithm to approximate MLEs obtained from the NONMEM package.

Simulated data are given by:

y ij = φ 1i - φ 2i x j φ 3i + x j + ε ij ; 1 ≤ i ≤ n , 1 ≤ j ≤ m (4.6)
where n = 30, m = 6, x 1 = 0, x 2 = 5, x 3 = 10, x 4 = 20, x 5 = 40 and x 6 = 80. The random effects and the additive noise are simulated with Gaussian distributions:

φ 1i ∼ iid N (105, 64) , φ 2i ∼ iid N (12, 36) , φ 3i ∼ iid N (10, 12.25) , ε ij ∼ iid N (0, 4).
According to [START_REF] Sheiner | A simulation study comparing designs for dose ranging[END_REF] and [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF], this model can be used for the analysis of blood pressure y as a function of the dose d of an anti-hypertensive drug from a longitudinal study. [START_REF] Walker | An EM algorithm for nonlinear random effects models[END_REF] compares different popular methods of estimation,such as FOCE (First-Order Conditional Estimation) and LAPLACIAN methods of NONMEM. He computes the means and the standard errors based on 50 simulations for these different estimators. Table 3 reproduces these values, with the estimates and standard errors obtained with the SAEM algorithm as well. We see that the EM algorithm of Walker and the SAEM algorithm give similar results, but it is important here to observe that only 300 iterations of the SAEM algorithm are performed for each single simulated dataset. Computing time for the SAEM algorithm is then very much reduced, in comparison to the Monte-Carlo EM algorithm that requires a sampling of 10,000 random variates at each iteration and converges very slowly. Table 3 also gives the estimation of the standard deviation of the MLE, using the approach proposed in Section 2.3. This method seems to be very accurate since these values are just below the empirical standard deviation computed from the 50 simulations. It should also be recalled that the gap observed between the true value of the parameter and the limit of the parameter estimate sequence generated by the algorithm results from two phenomena: first, the convergence of the maximum likelihood estimate toward the true value of the parameter that occurs as the number of individuals n tends to infinity; second, the convergence of the parameter estimate sequence generated by the algorithm toward a local maximum of the observed likelihood that occurs as the number of iterations of the algorithm tends to infinity.

Statistical modeling and estimation for Deformable Template Models

In this section, I present my contributions to the analysis of deformable template models, from modeling to parameter estimation.

The Bayesian Mixed Effects Template Model

We considered the hierarchical Bayesian framework for dense deformable templates developed in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] . Each image of a given population is assumed to be generated as a noisy and randomly deformed version of a common template drawn from a prior distribution on the set of templates. Individual deformations are hidden variables of the model or equivalently random effects in the mixed effects setting, whereas the template and the law of the deformations are parameters (or equivalently fixed effects) of interest.

The observation model

We considered gray level images (y i ) 1≤i≤n observed on a grid of pixels

{v u ∈ D ⊂ R 2 , u ∈ Λ} that is embedded in a continuous domain D ⊂ R 2 , (typically D = [-1, 1] × [-1, 1]
). Although the images are observed only at the pixels (v u ) u , we are looking for a template image I 0 : R 2 → R defined on the plane (the extension to images on R d is straightforward). Each observation y is assumed to be the discretization on a fixed pixel grid of a deformation of the template, plus independent noise. For each observation, there exists an unobserved deformation field z : R 2 → R 2 such that for u ∈ Λ:

y(u) = I 0 (v u -z(v u )) + (u) ,
where denotes an independent additive noise.

Considering the template and the deformations as continuous functions would lead to a dense problem.

We use the same framework as chosen in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] to describe both the templates I 0 and the deformation fields z. Our model takes two complementary aspects into account: photometric -indexed by p-corresponding to the templates and the noise variances, and geometric -indexed by g-corresponding to the deformations. We choose a representation of both the templates I 0 and the deformations z by finite linear combinations of the kernels centered at some fixed landmark points in the domain D: (v p,j ) 1≤j≤kp (respectively, (v g,j ) 1≤j≤kg ). They are therefore parameterized by the coefficients α ∈ R kp and β ∈ (R kg ) 2 that yield: ∀v ∈ D,

I α (v) (K p α)(v) kp j=1 K p (v, v p,j )α j , z β (v) (K g β)(v) kg j=1 K g (v, v g,j )β j .
For the sake of clarity, we denote the collection of data and their corresponding deformation coefficients by y t = (y t 1 , . . . , y t n ) and β t = (β t 1 , . . . , β t n ), respectively. The statistical model of the observations we consider is a generative hierarchical one. We assume conditional normal distributions for y and β:

     β ∼ ⊗ n i=1 N 2kg (0, Γ g ) | Γ g , y ∼ ⊗ n i=1 N |Λ| (z βi I α , σ 2 Id) | β, α, σ 2 , (5.1)
where ⊗ denotes the product of distributions of independent variables and zI α (u) = I α (v u -z(v u )), for u in Λ denotes the action of the deformation on the template image. The parameters of interest are α which determines the template image; σ 2 , the variance of the additive noise; and Γ g , the covariance matrix of the variables β. We assume that θ = (α, σ 2 , Γ g ) belongs to an open parameter space Θ:

Θ { θ = (α, σ 2 , Γ g ) | α ∈ R kp , α < R |, σ > 0, Γ g ∈ Sym + 2kg } ,
where . is the Euclidean norm, Sym + 2kg is the cone of real positive 2k g ×2k g definite symmetric matrices, and R is an arbitrary positive constant.

The Bayesian Statistical Model

Even though the parameters are finite dimensional, the maximum likelihood estimator can yield degenerate estimates when the training sample is small. Introducing prior distributions on the parameters regularized the estimation with small samples. The analytical effect of such priors can be seen in the parameter update steps (cf. [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]). We use a generative model based on standard conjugate prior distributions for parameters θ = (α, σ 2 , Γ g ) with fixed hyper-parameters. Specifically, we assume a normal prior for α, an inverse-Wishart prior on σ 2 and an inverse-Wishart prior on Γ g . Furthermore, all priors are assumed to be independent. This yields θ = (α, σ 2 , Γ g ) ∼ q para ν p ⊗ ν g where

         ν p (dα, dσ 2 ) ∝ exp - 1 2 (α -µ p ) t (Σ p ) -1 (α -µ p ) exp - σ 2 0 2σ 2 1 √ σ 2 ap dσ 2 dα, a p ≥ 3 , ν g (dΓ g ) ∝ exp(-Γ -1 g , Σ g F /2) 1 |Γ g | ag dΓ g , a g ≥ 4k g + 1 .
(5.2)

For two matrices A and B, we define A, B F tr(A t B) as the Frobenius dot product on the set of matrices where tr denotes the trace of the matrix.

Maximum A Posteriori Estimation

In this Bayesian framework, the parameter estimation will be performed by Maximum A Posteriori (MAP) defined as:

θn = argmax θ∈Θ p(θ|y) ,
where p denotes the posterior likelihood of the parameters given the observations. The dependence on n refers to the sample size. The existence and consistency (as the number of observed images tends to infinity) has been proven (see [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]). This contrasts with earlier studies in [START_REF] Glasbey | A penalised likelihood approach to image warping[END_REF] because it uses a penalized likelihood or the more recent maximum description length approach in [START_REF] Marsland | A minimum description length objective function for groupwise non rigid image registration[END_REF] for which consistency cannot be proved because the deformations are considered as nuisance parameters to be estimated.

However, the maximization problem of the posterior distribution has no closed form in our case, which prevents a direct computation of θn . The EM algorithm, although quite natural for maximizing a likelihood under a hierarchical model with missing variables, is not adapted to the deformable template model. In fact, direct computation is unfortunately not tractable and we have to find a solution to overcome the problematic E step where we have to compute an expectation with respect to the conditional distribution of β given y.

A first attempt was proposed in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] where this conditional distribution is approximated by a Dirac distribution at its mode. The authors called their algorithm the Fast Approximation with Mode EM (FAM-EM). The results were very interesting, but the authors point out the lack of convergence of the FAM-EM algorithm when the quality of the input images is not good, and, typically, when they are noisy. Note that the FAM-EM algorithm mentioned in that paper corresponds analytically exactly to the EM-Laplace (see [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF]). Indeed, the algorithm achieves its limits in several cases. Using the SAEM algorithm would make it necessary to sample the hidden variable from the conditional distribution. This sampling is not possible in this complex model.

We therefore apply the SAEM-MCMC algorithm with truncation on random boundaries proposed in [A3] to approximate the MAP estimator θn .

But the setting of deformable template models deals with high-dimensional missing variables. This raises several issues. If we simulate candidates for the hidden variable as a complete vector, it appears that most of the candidates are rejected. This is a typical high dimensional concentration phenomenon: locally, around a current point, the proportion of the space occupied by acceptable moves becomes negligible when the space dimension grows. From a more practical point of view, even if the proposed candidate is drawn with respect to the current prior distribution, it creates a deformation that is very different from the current one and too large for the corresponding deformed template to fit the observations. This yields very few possible moves from the current missing variable value, and the algorithm is stuck in a non-optimal location or converges very slowly. As a consequence, the infinite support of the unobserved variables is not well covered by the simulated variables.

One solution is to update the chain one coordinate at a time, conditional on the others. This corresponds to a Gibbs sampler and leads to more relevant candidates that have a higher chance of being accepted (see [START_REF] Amit | Convergence properties of the Gibbs sampler for perturbations of Gaussians[END_REF]). From an image analysis point of view, this puts stronger conditions on the type of deformations that are produced when proposing a candidate for each coordinate. Knowing the tendency of the movement given by the other coordinates, the candidate will either confirm it or not depending on if this is a suitable movement. It will thus be accepted with a corresponding probability. Even if some coordinates remain unchanged, some others are updated, enabling the algorithm to visit a larger part of the missing variable support. However, as the dimension grows, the computation time becomes very long. We thus adopt another approach: instead of updating only one coordinate, we optimize the proposal distribution in a high dimension. This can be efficiently done by using the AMALA sampler presented in Sections 3.1 and 3.2.

Estimation in the Bayesian mixture model

In collaboration with Stéphanie Allassonnière, we considered the framework of the multicomponent model introduced in Allassonnière et al. [2007] and extended the estimation method developed for the deformable template model to this multicomponent model [A5].

Let us consider the statistical estimation of the component weights and of the image labels in a multicomponent model given a set of images. In the existing methods, the templates of each component and the label are estimated iteratively (for example in methods like K-means) but the geometry, and related to this the metric used to compute the distances between elements, is still fixed. Moreover, the label, which is not observed, is, as the deformations, considered as a parameter and not as a hidden random variable. These methods do not lead to a statistical coherent framework for the understanding of deformable template estimation and none of these iterative algorithms derived from those approaches have a statistical interpretation in terms of parameter optimisation in a generative model that describes the data.

We therefore considered the statistical framework for dense deformable templates developed in Allassonnière et al. [2007] in the generalized case of a mixture model for multicomponent estimation. Each image of a database is assumed to be generated from a noisy random deformation of a template image picked randomly among a given set of possible templates. All the templates are assumed to be drawn from a common prior distribution on the template image space. To propose a generative model, each deformation and each image label have to be considered as hidden variables. The templates, the parameters of the deformation laws and the components weights are the parameters of interest. This generative model allows us to automatically break down the database into components and, at the same time, to estimate the parameters corresponding to each component while increasing the likelihood of the observations.

We consider a mixture of the deformable template models, which allows a fixed number τ m of components in each training set. This means that the data will be separated into τ m (at most) different components by the algorithm. For each observation y i , we consider the pair (β i , τ i ) of unobserved variables that correspond to the deformation field and to the label of image i, respectively. We denote this below by y t (y t 1 , . . . , y t n ), by β t (β t 1 , . . . , β t n ) and by τ t (τ 1 , . . . , τ n ). The generative model is:

                 τ ∼ ⊗ n i=1 τm t=1 ρ t δ t | (ρ t ) 1≤t≤τm , β ∼ ⊗ n i=1 N (0, Γ g,τi ) | τ , (Γ g,t ) 1≤t≤τm , y ∼ ⊗ n i=1 N (z βi I ατ i , σ 2 τi Id |Λ| ) | β, τ , (α t , σ 2 t ) 1≤t≤τm , (5.3)
where

z β I α (u) = I α (v u -z β (v u ))
is the action of the deformation on the template I α , for u in Λ, and δ t is the Dirac function on t. The parameters of interest are the vectors (α t ) 1≤t≤τm coding the templates, the variances (σ 2 t ) 1≤t≤τm of the additive noises, the covariance matrices (Γ g,t ) 1≤t≤τm of the deformation fields and the component weights (ρ t ) 1≤t≤τm . We denote the parameters by (θ t , ρ t ) 1≤t≤τm so that θ t corresponds to the parameters composed of the photometric part (α t , σ 2 t ) and the geometric part Γ g,t for component t. We assume that for all 1 ≤ t ≤ τ m , the parameter θ t = (α t , σ 2 t , Γ g,t ) belongs to the open space Θ defined as Θ

= { (α, σ 2 , Γ g ) | α ∈ R kp , |α| < R, σ > 0, Γ g ∈ Sym +
2kg } , where R is an arbitrary positive constant and Sym + 2kg is the set of strictly positive symmetric matrices. Concerning the weights (ρ t ) 1≤t≤τm , we assume that they belong to the set = (ρ

t ) 1≤t≤τm ∈]0, 1[ τm | τm t=1 ρ t = 1 .
We choose a normal distribution for the unobserved deformation variable because of the background we have in image analysis. In fact, the registration problem is an issue that has been studied in depth over the past two decades. The goal is, given two images, to find the best deformation that will match one image close to the other. Such methods require choosing the kind of deformations that are allowed (smooth, diffeomorphic, etc). These conditions are equivalent, for some of these methods, to choose a covariance matrix that enables to define an inner product between two deformations coded by a vector β (cf. [START_REF] Amit | Structural image restoration through deformable templates[END_REF], [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF]). The regularization term of the matching energy in the small deformation framework treated in this paper can be written as: β t Γ -1 g β. This looks like the logarithm of the density of a Gaussian distribution on β with 0 mean and a covariance matrix Γ g . The link between these two points of view has been given in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]; the mode of the posterior distribution equals the solution of a general matching problem. This is why we set such a distribution on the deformation vector β. Moreover, many experiments have been run using a large variety of such a matrix that now gives us a good initial guess for our parameter. This leads us to consider a Bayesian approach with a weakly informative prior.

We use a generative model that includes natural standard conjugate prior distributions with fixed hyper-parameters. These distributions are an inverse-Wishart prior on each Γ g,t and σ 2 t and a normal prior on each α t , for all 1 ≤ t ≤ τ m . All priors are assumed to be independent.

For the prior law ν ρ , we choose the Dirichlet distribution, D(a ρ ), with density:

ν ρ (ρ) ∝ τm t=1 ρ t aρ
, with fixed parameter a ρ .

For the sake of simplicity, let us denote this by N 2nk g and by T {1, . . . , τ m } n so that the missing deformation variables take their values in R N and the missing labels in T . We also introduce the following notations: η = (θ, ρ) with θ = (θ t ) 1≤t≤τm and ρ = (ρ t ) 1≤t≤τm .

In our Bayesian framework, we choose the MAP estimator to estimate the parameters: ηn = argmax η p(η|y) ,

(5.4)

where p(η|y) denotes the posterior distribution of η given the observations y.

In practice, to reach this estimator, we maximize this posterior distribution using a Stochastic Approximation EM (SAEM) algorithm coupled with a Monte Carlo Markov Chain (MCMC) method as used for the one component model in [A3]. However in the multicomponent model we show that it cannot be driven numerically. In fact, the direct generalization of the algorithm presented in [A3] turns out to be of no use in practice because of some trapping state problems.

If we consider the full vector (β, τ ) as a single vector of missing data, we can use the hybrid Gibbs sampler on R N × T . Even if this procedure provides an estimated parameter sequence which would theoretically converge toward the MAP estimator, in practice, as mentioned in [START_REF] Robert | Méthodes de Monte Carlo par chaînes de Markov[END_REF], it would take a quite long time to reach its limit because of the trapping state problem: when a small number of observations are assigned to a component, the estimation of the component parameters is hardly concentrated and the probability of changing the label of an image to this component or from this component to another is really small (most of the time under the computer precision).

We can interpret this from an image analysis viewpoint: the first iteration of the algorithm gives a random label to the training set and computes the corresponding maximiser η = (θ, ρ). Then, for each image, according to its current label, it simulates a deformation field which only takes into account the parameters of this given component. Indeed, the simulation of β through the Gibbs sampler involves a proposal whose corresponding Markov chain has q(β|τ , y, η) as stationary distribution. Therefore, the deformation tries to match y to the deformed template of the given component τ . The deformation field tries to get a better connection between the component parameters and the observation, and there is only small probability that the observation given this deformation field will be closer to another component. The update of the label τ is therefore conditional to this deformation which would not leave much chance to switch component.

To overcome the trapping state problem, we will simulate the optimal label, using as many Markov chains in β as the number of components so that each component has a corresponding deformation which "computes" its distance to the observation. Then we can simulate the optimal deformation corresponding to that optimal label.

Since we aim to simulate (β, τ ) through a transition kernel that has q post (β, τ |y, η) as stationary distribution, we simulate τ with a kernel whose stationary distribution is q post (τ |y, η) and then β through a transition kernel that has q post (β|τ , y, η) as stationary distribution.

For the first step, we need to compute the weights q(t|y i , η) ∝ q(t, y i |η) for all 1 ≤ t ≤ τ m and all 1 ≤ i ≤ n which cannot be easily reached. However, for any density function f , for any image y i and for any 1 ≤ t ≤ τ m , we have

q(t, y i |η) = E qpost(β|yi,t,η) f (β) q(y i , β, t|η) -1
.

(5.5)

Obviously the computation of this expectation w.r.t. the posterior distribution is not tractable either but we can approximate it by a Monte Carlo sum. However, we cannot easily simulate variables through the posterior distribution q post (•|y i , t, η) as well, so we use some realisations of an ergodic Markov chain having q post (•|y i , t, η) as stationary distribution instead of some independent realisations of this distribution.

The solution we propose is the following: suppose we are at the k th iteration of the algorithm and let η be the current parameters. Given any initial deformation field ξ 0 ∈ R 2kg , we run, for each component t, the hybrid Gibbs sampler Π η,t on R 2kg J times so that we get J elements ξ t,i = (ξ (l) t,i ) 1≤l≤J of an ergodic homogeneous Markov chain whose stationary distribution is q(•|y i , t, η). Let us denote by ξ i = (ξ t,i ) 1≤t≤τm the matrix of all the auxiliary variables. We then use these elements for the computation of the weights p J (t|ξ i , y i , η) through a Monte Carlo sum:

p J (t|ξ i , y i , η) ∝ 1 J J l=1 f (ξ (l) t,i ) q(y i , ξ (l) t,i , t|η) -1 , (5.6)
where the normalisation is done such that their sum over t equals one, involving the dependence on all the auxiliary variables ξ i . The ergodic theorem ensures the convergence of our approximation toward the expected value. We then simulate τ through ⊗ n i=1 τm t=1 p J (t|ξ i , y i , η)δ t .

Concerning the second step, we update β by re-running J times the hybrid Gibbs sampler Π η,τ on R N starting from a random initial point β 0 in a compact subset of R N . The size of J will depend on the iteration k of the SAEM algorithm in a sense that will be precised later, thus we now index it by k.

The density function f involved in the Monte Carlo sum above needs to be specified to get the convergence result proved in the last section of this paper. We show that using the prior on the deformation field enables to get the sufficient conditions for convergence. This density is the Gaussian density function and depends on the component we are working with:

f t (ξ) = 1 √ 2π 2kg |Γ g,t | exp - 1 2 ξ t Γ -1 g,t ξ .
(5.7)

Let (η k ) be the sequence generated by the algorithm 2.2 using the sampler described above. We obtained the following convergence result: Theorem 5.1. Under the assumptions of Theorem 4.2 of [A5], we have for all β 0 ∈ K, τ 0 ∈ T , s 0 ∈ S and η 0 ∈ Θ × , lim k→∞ d(η k , L) = 0 Pβ 0 ,τ 0 ,s0,0 -a.s , where Pβ 0 ,τ 0,s0 ,0, is the probability measure associated with the chain (Z k = (β k , τ k , s k , κ k )) k≥0 starting at (β 0 , τ 0 , s 0 , 0) and L { η ∈ η(S), ∂l ∂η (η) = 0}. Proof. Even if the only algorithmic difference between our algorithm and the SAEM algorithm is the simulation of the missing data which is not done with respect to the conditional law q post (β, τ |y, η) but through an approximation which can be arbitrarily close, this yields a very different proof. Indeed, whereas for the SAEM algorithm, the stochastic approximation leads to a Robbins-Monro type equation with no residual term r k , our method induces one. The first difficulty is therefore to prove that this residual term tends to 0 while the number of iterations k tends to infinity. Our proof is decomposed into two part, the first one concerning the deformation variable β and the second one the label τ . The first term requires to prove the geometric ergodicity of the Markov chain in β generated through our kernel. For this purpose, we prove some typical sufficient conditions which include the existence of a small set for the transition kernel and a drift condition. Then, we use for the second term some concentration inequalities for non stationary Markov chains to prove that the kernel associated with the label distribution converges toward the conditional distribution q post (τ |y, η).

The second difficulty is to prove the convergence of the excitation term e k . This can be carried out as in [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF] using the properties of our Markov chain and some martingale limits properties.

Estimation in the Large Diffeomorphic Deformable Mapping Model

In collaboration with Stéphanie Allassonnière and Stanley Durrlemann, we applied the Anisotropic Metropolis Adjusted Langevin Algorithm as a sampler into the SAEM-MCMC algorithm to estimate the parameters in the Large Diffeomorphic Deformable Mapping Model [A13]. We also proposed an extension of the model to optimize the position of the control points as well as a criterion to select the number of control points.

Introduction

Formerly we consider a simple deformation model may be the so-called "linearized deformation". A linearized deformation φ is defined by the displacement field v of each point in the domain D ⊂ R d : ∀r ∈ D, φ(r) = r + v(r). The main advantage of this class of deformations is its numerical simplicity as it parameterizes the deformation by a single vector field v. Nevertheless, even with regularity conditions on v, there is no guarantee that the deformation is invertible, meaning that the deformation may create holes or overlapping regions in the domain. To avoid such unrealistic behaviors, diffeomorphic maps that preserve the topology of the shapes in the image set should be considered. This means that we assume that every sample has the same topology or, equivalently, that differences within sample shapes do not rely on changes in topology.

Diffeomorphic deformations can be built on linearized deformations within the framework of the Large Diffeomorphic Deformation Metric Mapping (LDDMM), which was introduced by [START_REF] Trouvé | Diffeomorphism groups and pattern matching in image analysis[END_REF] and [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF] and further developed by [START_REF] Allassonnière | Geodesic shotting and diffeomorphic matching via textured meshes[END_REF], [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], Glaunès et al. [2003], [START_REF] Holm | Soliton dynamics in computational anatomy[END_REF], [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF], [START_REF] Miller | Group action, diffeomorphism and matching: a general framework[END_REF]. In this framework, the above linearized deformations are considered to be infinitesimal deformations, and the vector field v is seen as an instantaneous velocity field. The composition of such deformations creates a flow of diffeomorphisms, which can be written at the limit as the solution of a differential equation. The set of such diffeomorphisms can be equipped with a group structure and a right-invariant metric, providing regularity on the driving velocity fields. It follows that the set of images is given the structure of an infinite-dimensional manifold, on which distances are computed as the geodesic length in the deformation group between the identity map and the diffeomorphism that maps one image on to another.

It was shown by [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] that this infinite dimensional deformation set can be efficiently approximated by a finite control point parameterization that carries momentum vectors. This finite dimension reduction is a key aspect for statistical analysis. [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF] have enforced the velocity fields that are defined everywhere in the domain to be parameterized by a finite set of control points. Positions of control points are not given as a prior but optimized as parameters of the statistical model. As a consequence, control points tend to move to the regions that show the greatest variability among samples while optimizing a least-square criterion. At the same time, this optimization makes it possible to reduce the number of control points for the same matching accuracy, compared to the case where control points are fixed as the nodes of a regular lattice.

Once the deformation model has been fixed, it is necessary to estimate the parameters of the associated statistical model including, in particular, the template image. Different algorithms have been proposed to solve the template estimation. Most of them are based on a deterministic gradient descent. In particular, [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF] managed simultaneously the optimization in control point positions and momentum vectors using a joint gradient descent. Although it provided visually interesting results in several practical cases, the nature of the limit is not identified. Moreover, this type of method fails in specific cases, in particular, when using noisy training data. Another point of view is to consider stochastic algorithms, e.g., [START_REF] Zhang | Bayesian estimation of regularization and atlas building in diffeomorphic image registration[END_REF] using a Hamiltonian Monte Carlo sampler into a Monte Carlo Expectation Maximization algorithm in the dense LDDMM setting, although there is no theoretical convergence property proved for this algorithm.

We considered now the LDDMM setting where the deformations are parameterized by a finite number of initial control point positions and momenta such as in [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF]. To do this, we extend the generative statistical model of [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]. In that model, the deformations are assumed to be linearized and are modeled as random variables that are not observed. This enables us to estimate the representative parameters of their distribution that will characterize the geometric variability. On the one hand, we extend this approach to the LDDMM framework. On the other, we introduce the control point positions as population parameters into the model so that they can be optimized in the estimation process. This enables us to better fit the deformation model leading to a more accurate estimation of the geometric parameters.

From an algorithmic point of view, we propose to use the Anisotropic Metropolis Adjusted Langevin Algorithm (AMALA) within the SAEM algorithm introduced in [A12]. This algorithm has shown itself to have very interesting theoretical and numerical properties. Indeed, the AMALA sampler enables us to better explore the target distribution support in very high dimensions, compared to other samplers. It also increases the speed of convergence of the estimation algorithm. Moreover, we take advantage in our sampler of the efficient computation used in the joint gradient descent by [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF] so that the optimization of control point positions is of no additional cost at each iteration.

Another interesting question is how to optimize the number of control points required to parameterize the deformations. Indeed, the number of control points essentially depends on the variability in the data: it should be estimated rather than fixed by the user. In the geometric approach [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], control points were automatically selected using a L 1 type penalty that tends to zero out momentum vectors of small magnitude. Numerically it is solved by an adapted gradient descent known as FISTA (see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]). However, this penalty acts on each observation separately, meaning that a control point that is needed to match only a single observation will be kept in the final set of control points. From a statistical point of view, this control point can be thought of as an outlier that could preferably be removed from the basis. The L 1 penalty is also not suitable for statistical purposes, since its associated distribution, namely the Laplace prior, does not generate sparse variables. In other words, the criterion with the L 1 penalty that is minimized in [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF] could not be interpreted as the log likelihood of a statistical model that generates sparse solutions.

We propose to include a sparsity constraint in the parameter space of our statistical model through a thresholding step, borrowing ideas from the Group LASSO literature initiated in [START_REF] Bach | Consistency of the group lasso and multiple kernel learning[END_REF]. This has the advantage to select control points based on their importance for the description of the variability of the whole population, and not only of one single sample. The thresholding step is then included in the Maximization step, so that the same AMALA-SAEM algorithm can be used for the estimation process. We also exhibit a criterion to select an optimal threshold leading to an optimal number of control points.

The LDDMM model

The model of diffeomorphic deformations we choose is derived from the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework (see [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF]], Miller et al. [2002], [START_REF] Trouvé | Diffeomorphism groups and pattern matching in image analysis[END_REF]), which generalizes the linearized deformation setting that has been used for the statistical estimation of atlases by [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]. In the linearized deformation setting, the deformation φ is given by :

φ(r) = r + v(r), ∀ r ∈ D ,
(5.8) with d = 2 or 3, and v a vector field on R d .

To build a diffeomorphic map, we use the linearized deformations given in Equation (5.8) as infinitesimal steps, and consider the corresponding vector field as an instantaneous velocity field. More precisely, we consider time-dependent velocity fields (v t ) t for a time-parameter t varying in [0, 1]. The motion of a point r 0 in the domain of interest D describes a curve t → r(t) which is the integral curve of the following Ordinary Differential Equation (ODE) called flow equation :

   dr(t) dt = v t (r(t))
r(0) = r 0 .

(5.9)

The deformation φ 1 is defined as follows:

∀r 0 ∈ D, φ 1 (r 0 ) = r(1) .
Conditions under which this map φ 1 is diffeomorphic can be found in [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]. In particular, the existence, uniqueness and diffeomorphic property of the solution are satisfied if the velocity v t belongs to a RKHS at all time t and is square integrable in time.

Under these conditions, the model builds a flow of diffeomorphic deformations φ t : r 0 -→ r(t) for all t ∈ [0, 1]. The flow describes a curve in a sub-group of diffeomorphic deformations starting at the identity map. The RKHS V plays the role of the tangent space of such an infinite-dimensional Riemannian manifold at the identity map Id. We can provide this group of diffeomorphisms with a right-invariant metric, where the square distance between the identity map Id = φ 0 and the final deformation φ 1 is given as the total kinetic energy used along the path: d(Id, φ 1 ) 2 = 1 0 v t 2 V dt, where • V is the norm in the RKHS. The existence and uniqueness of minimizing paths have been shown by [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF].

According to mechanical principles, one can show that the kinetic energy is preserved along the geodesic paths, namely for all t ∈ [0, 1] v t V = v 0 V . Moreover, the velocity fields (v t ) along such paths satisfy Hamiltonian equations, meaning that the geodesic is fully parametrized by the initial velocity field v 0 . This velocity field plays the role of the Riemannian logarithm of the final diffeomorphism φ 1 . Therefore, it belongs to a vector space and allows the definition of tangent-space statistics in the spirit of [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF] and [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF].

Following [START_REF] Durrleman | Optimal data-driven sparse parameterization of diffeomorphisms for population analysis[END_REF] and [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], we further assume that v 0 is the interpolation of momentum vectors (α k,0 ) k at control point positions (v g,j ,0 ) k :

v 0 (r) = kg k=1 K g (r, v g,j ,0 )α k,0 ,
(5.10) where K g is the kernel associated to the RKHS V . In this context, it has been shown in [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF] that the velocity fields (v t ) t along the geodesic path starting at the identity map in the direction of v 0 keep the same form: (5.11) where the control point positions (v g,j (t)) k and the momentum vectors (α k (t)) k satisfy the Hamiltonian equations:

v t (r) = kg k=1 K g (r, v g,j (t))α k (t) ,
               dv g,j (t) dt = kg l=1 K g (v g,j (t), c l (t))α l (t) dα k (t) dt = -   kg l=1 d c k (t) (K g (v g,j (t), c l (t))α l (t))   t α k (t)
(5.12)

with initial conditions v g,j (0) = c 0,k and α k (0) = α 0,k for all 1 ≤ k ≤ k g . This is similar to the equations of motion of a set of k g self-interacting particles, with K g modeling the interactions. One can easily verify that the Hamiltonian defined as

H t = v t 2 V = kg k=1 kg l=1 α l (t) t K g (c l (t), v g,j (t))α k (t)
is constant in time when control point positions and momentum vectors satisfy the system (5.12). This model defines a finite dimensional subspace of the group of diffeomorphisms. For a given set of initial control points, the diffeomorphisms are parametrized by the momentum vectors attached to them. For one instance of the initial momentum vectors, one builds the motion of the control points and of the momentum vectors by integrating the Hamiltonian system (5.12). Then, they define a dense velocity field at each time t according to Equation (5.11). Finally, one can find the motion φ t (r 0 ) of any point r 0 in the domain D by integrating the flow equation (5.9). In this framework, the tangent-space representation of the diffeomorphic deformation φ 1 is given by the initial velocity field v 0 parametrized by z = ((c 0,k , α 0,k )) k , called the initial state of the particle system. The position φ 1 (r) depends on the parameters ((c 0,k , α 0,k )) k via the integration of two non-linear differential equations in Equation (5.12) and Equation (5.9).

Remark 5.1. The LDDMM framework formulation involves a coupling on the control point and the momentum evolutions along the geodesic path which is not the case in the linearized deformation setting. This joint equation introduces more constraints reducing the dimension of the solution space. Therefore, the identifiability of the control point positions may be expected in our LDDMM framework. This property would most probably fail in the linearized deformation setting where the momenta and the control points are not coupled.

Statistical model

As pointed out in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF], the gradient descent optimization with respect to the template together with the momenta does not necessarily converge if the training set is noisy. To solve this problem, we introduce here a statistical model where we consider the deformations as well as the control point positions as non-observed random variables, in the spirit of the BME template model [A3].

We choose to model our data by a generative hierarchical model. In this model, the distribution of the deformations in the diffeomorphism group is parametrized. In a statistical approach, these parameters are estimated from the data, thus providing a metric in the shape space which is adapted to the data and takes into account the deformation constraints. This is in contrast to geometric approaches that estimate the template using a fixed metric.

More precisely, let I 0 be a template image: I 0 : R d → R. We consider an observation, namely an image y, as a noisy discretization on a regular grid Λ of a diffeomorphic deformation of the template image. Let φ z 1 be the solution of both the flow equation (5.9) and the Hamiltonian system (5.12) with initial condition z = ((c 0,k , α 0,k )) k . Then, for all s ∈ Λ, (5.13) where σ denotes an additive centered Gaussian random noise on the grid Λ with variance σ 2 , and r s is the coordinate of the voxel s in the continuous domain D.

y(s) = I 0 ((φ z 1 ) -1 (r s )) + σ (s),
We are provided with n images y = (y i ) 1≤i≤n in a training set. We assume that each of them follows the probabilistic model (5.13) and that they are independent.

We consider the initial state of particles, namely the control point positions and the momentum vectors, as random variables and estimate their probabilistic distributions, restricting ourselves to the case of parametric distributions. We assume that control points live in the template domain D and that they are the same for all observations. By contrast, the momentum vectors attached to them are specific to each observation, as they parametrize the matching of the template with each sample image.

Therefore, we propose the following probabilistic model: we assume that the initial control point positions c 0 = (c 0,k ) 1≤k≤kg are drawn through a Gaussian distribution with mean c0 and covariance a c Id where Id is the identity matrix of dimension dk g . We define the initial momenta α 0 = (α i 0 ) 1≤i≤n with α i 0 = (α i 0,k ) 1≤k≤kg . We assume that the variables (α i 0 ) 1≤i≤n are independent identically distributed and follow a Gaussian distribution with mean 0 and covariance matrix Γ g . Note that this covariance matrix depends on the initial control point positions as the momenta are attached to them. Moreover the momenta α 0 are assumed to be independent of the control point positions c 0 given Γ g .

Following the same lines as [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF], we parametrize the template function I 0 as a linear combination of gray level values of fixed voxels (v p,j ) 1≤k≤kp equidistributed on the domain D. The interpolation kernel is denoted by K p and the combination weights are denoted by w: ∀r ∈ D,

I 0 (r) = kp k=1 K p (r, v p,j )w k .
(5.14)

The action of a diffeomorphism on this template is the linear combination of the deformed kernel with the same weights: ∀r ∈ D,

K z p w(r) = I 0 • (φ z 1 ) -1 (r) = kp k=1 K p (φ z 1 ) -1 (r), v p,j w k .
(5.15)

The parameters of the model are θ = (w, σ 2 , Γ g , c0 ) and the random variables (α 0 , c 0 ) are considered as hidden random variables. As we often deal with small sample size in practice, we restrict our inference to a Bayesian setting. Some of the priors can be informative as the one of Γ g . Other priors may be non-informative as for the expectation of the control point positions for which no additional information is available. The complete model writes therefore:

                         θ = (w, σ 2 , Γ g , c0 ) ∼ ν p ⊗ ν g c 0 ∼ N dkg (c 0 , a c Id)| θ, α i 0 ∼ N dkg (0, Γ g )| θ, ∀1 ≤ i ≤ n , y i ∼ N |Λ| (K (c0,α i 0 ) p w, σ 2 Id) | (c 0 , α i 0 ), θ, ∀1 ≤ i ≤ n .
(5.16)

We define the prior distributions as follows:

     ν g (dΓ g , dc 0 ) ∝ exp(-Γ -1 g , Σ g F /2) 1 √ | det(Γg)| ag • exp -1 2 (c 0 -µ c ) t Σ -1 c (c 0 -µ c ) dΓ g dc 0 , ν p (dw, dσ 2 ) ∝ exp -1 2 w t Σ -1 p w • exp - σ 2 0 2σ 2 1 √ σ 2 ap dwdσ 2 ,
where ., . F designs the Frobenius scalar product and the hyper-parameters satisfy a g ≥ 4k g +1, Σ g = Id, σ 2 0 > 0, a p ≥ 3 and Σ p is derived from the interpolation kernel K p and the photometric grid (v p,j ) 1≤k≤kp (see [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] for more details). Concerning the hyper-parameters of the control point prior (µ c , Σ c ), we choose µ c to be the vector of the equidistributed grid coordinates. The covariance matrix Σ c is assumed non-informative. All priors are the natural conjugate priors and are assumed independent to ease derivations.

Remark 5.2. From a modeling point of view, the positions of the control points could have been consider as parameters of our model since they are fixed effects of the whole population as well as the template. However considering control points as parameters does not lead to a model belonging to the exponential family. Thus, we could not benefit from the convergence properties and efficient implementation of the SAEM algorithm for this class of models. Therefore, we model the control point positions as random variables following a Gaussian distribution.

Parameter estimation

Let us define y = (y 1 , ..., y n ). We consider the Maximum A Posteriori (MAP) estimator denoted by θn obtained by maximizing the posterior density of θ conditional to y as follows : θn = argmax θ p(θ|y).

(5.17)

We first show that for any finite sample the maximum a posteriori will lie in the parameter set Θ; this is non-trivial due to the highly non-linear relationship between parameters and observations in the model.

Theorem 5.2 (Existence of the MAP estimator). For any sample y, there exists θn ∈ Θ such that q( θn |y) = sup θ∈Θ q(θ|y) .

We are interested in the consistency properties of the MAP estimator without making strong assumptions on the distribution of the observations y denoted by P . We seek to prove the convergence of the MAP estimator to the set Θ * defined by :

Θ * = { θ * ∈ Θ | E P (log q(y|θ * )) = sup θ∈Θ E P (log q(y|θ))}.
Theorem 5.3 (Consistency). Assume that Θ * is non empty. Then, for any compact set K ⊂ Θ, for all ε > 0, lim

n→+∞ P ( δ( θn , Θ * ) ≥ ε ∧ θn ∈ K ) = 0 ,
where δ is any metric compatible with the usual topology on Θ.

The proof follows the lines of [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]. Indeed, the observed likelihood of our diffeomorphic BME template model has the same regularity properties and asymptotic behaviors in the parameters as the linearized one.

To compute the MAP, we use the AMALA-SAEM algorithm proposed in [A12]. Indeed in our applications, the missing variables composed of the initial momenta and positions of control points z = (c 0 , α 0 ) are of very high dimension.

Extension toward sparse representation of the geometric variability

Obviously, the number of degrees of freedom needed to describe the variability of a given shape should be adapted to this shape. Therefore, the number of control points in our model should be estimated as a parameter of the model and not fixed by the user. This leads to automatically optimize the dimension of the deformation model. We propose here to simultaneously optimize the positions of the control points and select a subset of the most relevant ones for the description of the variability.

In [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], the control point selection is done adding an L 1 penalty on the momenta to the energy E θ and performing an adapted gradient descent called FISTA (see [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]). The effect of this penalty is to zero out momenta of small magnitude and to slightly decrease the magnitude of the other ones. A control point which does not contribute to at least one of the templateto-observation deformations at the convergence of the algorithm is called inactive. Note that since control points move in the domain, inactive control points may become active during the optimization process, and vice-versa.

This method suffers from three main limitations. First, the Laplace prior associated to the L 1 penalty does not generate sparse observations. Second, the method keeps active control points that may contribute to only few template-to-observation deformations. Lastly, L 1 penalty implies a soft thresholding step on the momentum vectors, thus reducing the norm of these vectors keeping the direction and therefore the local curvature. As a consequence, important momenta for the description of the variability will also be penalized. In the following, we propose to select control points given their importance to describe the variability of the whole population, and not of outliers. The idea is to inactivate a control point if the distribution of the momenta attached to it is not strongly correlated with the momentum distribution of other control points. Therefore our procedure selects control point positions and their number, relevant with regards to the whole population. This constraint on the momenta is taken into account in the model by assuming that the geometric covariance matrix Γ g is of the form Γ g = A g + ε g Id, where ε g is a small positive real number and A g is a sparse symmetric positive matrix. We introduced a positive threshold λ which plays an equivalent role as the weight of the L 1 penalty in the criterion optimized in [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF]. The larger, the sparser the solution. We do not consider a control point whose contribution to A g is lower than the threshold in a given sense defined through the parameter update equation [A13] This modified update is performed at each iteration of the estimation algorithm in the M-step.

To go one step further, we propose to automatically select an optimal threshold λ. We consider a criterion based on two relevant quantities namely the data attachment residual over the training images and the number of active control points. Indeed, the larger the threshold, the larger the residual and the lower the number of active control points. These quantities are computed for different values of the threshold. These sequences are then normalized to 1. The optimal threshold is chosen to be the point where the two normalized sequences intersect.

Experiments

We apply all the proposed estimation algorithms in the BME template models on different datasets. The first one is the USPS hand-written digit base as used in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF]. The other two are medical images of 2D corpus callosum [A4] and 3D murine dendrite spine excrescences [A12]. We present in this manuscript only the results obtained for the template estimation. The corresponding results for the geometrical parameter and the residual variance are given in the corresponding papers.

We begin with presenting the experiments on the USPS database. In order to make comparison, we estimate the parameters in the same conditions as in the previous mentioned works that is to say using the same 20 images per digit. Each image has grey level between 0 (background) and 2 (bright white). These images are presented on the left panel of Figure 2. We also use a noisy training dataset generated by adding a standardized independent Gaussian noise. These images are presented on the right panel of Figure 2. We test five algorithms: the deterministic approximation of the EM algorithm (FAM-EM) presented in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF], four SAEM-MCMC where the sampler is either the MALA, the adaptive MALA proposed in [START_REF] Atchadé | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF], the hybrid Gibbs sampler presented in [A3] and our AMALA algorithm.

Computational performances

The computational time is smaller for the three MCMC-SAEM algorithms using "MALA-like" samplers compared to the FAM. Comparing to the hybrid Gibbs-SAEM, the computational time is 8 times lower with the AMALA-SAEM in this particular case of application. Indeed, the hybrid Gibbs sampler requires no computation of the gradient. However, it includes a loop over the coordinates of the hidden variable, here the deformation vector of size 2k g = 72. At each of these iterations, the candidate is straightforward to sample whereas the computational cost lies into the acceptance rate. When this becomes heavy, the less times you calculate it, the better. In the AMALA-SAEM, this acceptance rate only has to be calculated once for each image. Therefore, even when the dimension of the hidden variable increases, this is of constant cost. The main price to pay is the computation of the gradient. Therefore, a tradeoff has to be found between the computation of either one gradient or dk g acceptance rates in order to select the algorithm to use.

Results on the template estimation

All the estimated templates obtained with the five algorithms and noise-free and noisy training data are presented in Figure 3. As noticed in [A3], the FAM-EM estimation is sharp when the training set is noise-free and is deteriorated while adding noise. This behavior is not surprising with regard to the theoretical bound established in Bigot and Charlier [2011] in the particular case of compact deformation group. Considering the adaptive sampler, it does not reach a good estimation of the templates which are still very blurry and noisy in both cases. The problem seems to come from the very low acceptation rate already at the beginning of the estimation. The bad initial guess we have about the covariance matrix of the proposal seems to block the chain. Moreover, the tuning parameters are difficult to calibrate along the iterations of the estimation algorithm. Concerning the estimated templates using the Gibbs, MALA and AMALA samplers, they look very similar to each other using the noise-free data as well as the noisy ones. This similarity confirms the convergence of all these algorithms toward the MAP estimator. In this case, the templates are as expected: noise free and sharp.

Nevertheless, when the dimension of the hidden variable increases, both the Gibbs and the MALA samplers show limitations. We run the estimation on the same noisy USPS database, increasing the number k g of geometrical control points. We choose the dimension of the deformation vector equal to 72, 128 and 200. The Gibbs-SAEM would produce sharp estimations but explodes the computational time. For this reason, we did not run this algorithm on higher dimension experiments. The results are presented in Figure 4. Concerning the MALA sampler, it does not seem to capture the whole variability of the population in such high dimension. This yields a poorly estimation of the templates. This phenomenon does not appear using our AMALA-SAEM algorithm. The templates still look sharp and the acceptation rates remain reasonable.

2D medical image template estimation

A second database is used to illustrate our algorithm. As before, in order to make comparisons with existing algorithms, we use the same database presented in [A4]. It consists of 47 medical images, each of them is a 2D square zone around the end point of the corpus callosum. This box contains a part of this corpus callosum as well as a part of the cerebellum. Ten exemplars are presented in the top rows of Figure 5.

The estimations are compared with these obtained with the FAM-EM and the hybrid Gibbs-SAEM algorithms and with the grey level mean image (bottom row of Figure 5). In this real situation, the Euclidean grey level mean image (a) is very blurry. The estimated template using the FAM-EM (b) provides a first amelioration in particular leading to a sharper corpus callosum. However, the cerebellum still looks blurry in particular when comparing it to the shape which appears in the template estimated using the hybrid Gibbs SAEM (c). The result of our AMALA-SAEM is given in image (d). This template is very close to (c) as we could expect at a convergence point. Nevertheless the AMALA-SAEM has much lower computational time than the hybrid Gibbs-SAEM. This shows the advantage of using AMALA-SAEM in real cases of high dimension.

Experiments in the LDDMM model

We consider the model with random control points presented in Section 5.4.3 as well as its simplified version where the control points are fixed. The number of control points is chosen equal to 4, 9 or 16 depending on the experiments. We infer the atlas of each digit independently using our stochastic estimation algorithm for the two models. We present the estimated templates obtained with both models and varying number of control points in Figure 6. The first row shows the template images estimated with control points fixed. The second one provides the estimated templates together with the estimated control point positions.

As expected, the contours in the template image become sharper in both cases as the number of control points is increased. Moreover, the number of control points being fixed, the sharpness of the estimated template is improved by allowing the control points to move toward optimized positions. We can also note that the estimated control points are informative as they tend to move toward the contours of the digits, and in particular toward those that correspond to the regions of highest variability among samples. It is particularly noticeable on digits 5 and 6 for example.

Mouse mandible experiment in LDDMM model

We consider a second training set composed of 36 X-ray scans of mouse mandibles. Five of them are presented in Figure 7. The estimated template images resulting from three different experiments are shown in Figure 8. The image on the left shows the template estimated using 260 fixed equidistributed control points. The image on the middle (resp. right) shows the estimated template using 117 (resp. 70) estimated control points. These templates look similar, thus showing that the same photometric invariants have been captured in each experiment. These invariants include the main bones of the mandibles (i.e. the brightest areas in the image). The decrease in number of control points is balanced by the optimization of their optimal positions. Control points in the right image are noticeably located on the edges of the shape in order to drive the dilation, contraction and opening of the mandible. Depending on the desired precision of the atlas, we can reduce even more the number of control points. This enables a faster estimation task at the cost of providing less information about the data.

Us Postal dataset experiment in LDDMM model

We conduct different experiments with different thresholds λ between 0.3 and 0.8 in order to see the evolution of the sparsity with respect to this parameter and also to capture the most interesting one (depending on the training digit). The initial number of control points is set to 16. The results of these experiments are presented in Figure 9. As expected, increasing the threshold λ decreases the final number of selected control points, whose effects on template sharpness and description of variability have been presented in Figure 6. Using the modified parameter update equation to enforce sparsity allows to automatically select a subset of control points leading to estimation results of the same accuracy (see Figure 9). Contrary to the L 1 prior used by [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF], our sparsity prior selects a small number of control points without penalizing the magnitude of the momenta. Hence the variability of the model is not under-estimated. In this respect, our thresholding process has an effect which is closer to the expected L 0 norm than its surrogate L 1 norm. Independently of the threshold λ, control points move in areas where the shape is the most variable. This can be noticed in the loop of the digit 2 which is highly variable, especially in contrast to the loop of the digit 6 which is much more stable in shape across observations. This can be seen by a fastest decrease in number of control points when the threshold λ is increased for the digit 6 compared to digit 2. It is also interesting to notice how our model deals with a mixture of 2 that could be written with or without a loop. Such variability violates the hypothesis of our model, which assumes that observations derive from a diffeomorphic deformations of the template image. In this situation, the model estimates a template image that is fuzzy in the region of the loop: the non-diffeomorphic variability has been interpreted as a photometric variation. To overcome this problem, one may investigate the use of several template images in the atlas along the lines of [A5]. The optimal threshold is chosen applying the proposed criterion. Figure 10 shows the estimated templates with their control points corresponding to the optimal threshold. The number of control points reflects the variability of the digits. In particular, very constrained shapes (see digits 1 and 9) require fewer control points than very complex irregular forms (see digits 3 and 8).

Experiments in the mixture model

We considered the multicomponent model decribed in Section 5.3. We ran the SAEM-MCMC algorithm presented in Section 5.3 on the US Postal dataset. In Figure 11, we showed the two estimated templates obtained in the mixture model by the algorithm with 40 training examples per class. It appeared that the two components reached were meaningful, such as the 2 with and without loop or American and European 7.

Goodness-of-fit test for Gaussian regression with block correlated errors

In collaboration with Sylvie Huet, we proposed a goodness-of-fit test for testing linear hypothesis on the expectation of a Gaussian vector with block correlated errors [A10].

Introduction

We consider n independent Gaussian vectors Y i , i = 1, . . . , n with unknown expectation f i and covariance matrix Σ i known up to some unknown parameters. The covariance matrix W of the random variable Y = (Y T 1 , . . . , Y T n ) T has a block diagonal structure composed of n squared blocks. This assumption means that the n blocks behave independently but that a correlation exists among the observations within the same block. Many models in various application fields may be covered by this framework, e.g., Gaussian multiresponse regression models with heteroscedastic errors, autocorrelated error models and mixed-effects models.

For the sake of simplicity, we assume that all Y i have the same size J. Let us denote the expectation of

Y by f = (f T 1 , . . . , f T n ) T .
Our aim is to test the null hypothesis that f belongs to a specific linear subspace V of R nJ against the alternative that it does not, without making any other assumption about its covariance matrix W than the block diagonal structure detailed above. For example, it would be possible to test the linearity of the relationship between the response and one covariate, or to test that a subset of covariates suffices to model the response.

Checking the adequacy of the model is of particular interest when analyzing data. Several graphical tools are proposed in regression models and, in particular, in linear mixed-effects models (see [START_REF] Diggle | Analysis of longitudinal data[END_REF], for example). Our interest is in goodness-of-fit testing procedures, also referred to as lack-of-fit or specification tests. Our objective is to propose a procedure without any a priori knowledge about f , which is adaptive rate-optimal and consistent against local alternatives.

A recent review on goodness-of-fit tests for regression models is provided by Gonzalez-Manteiga and Crujeiras [2013]. Several papers have considered the case where f i takes the particular form f i = (µ(X i,1 ), . . . , µ(X i,J )) T for each i = 1, . . . , n, for some function µ and for some covariates (X i,j , j = 1, . . . , J), taking values in R r . Some tests are based on nonparametric function estimation methods since it appears logical to compare a nonparametric estimation of µ to a parametric one computed under the null hypothesis. We refer, for example, to [START_REF] Chen | Testing for no effect in nonparametric regression via spline smoothing techniques[END_REF], [START_REF] Eubank | Testing for no effect in nonparametric regression[END_REF], [START_REF] Härdle | Comparing nonparametric versus parametric regression fits[END_REF], [START_REF] Härdle | Testing parametric versus semiparametric modeling in generalized linear models[END_REF], [START_REF] Hart | Nonparametric smoothing and lack-of-fit tests[END_REF], [START_REF] Müller | Goodness-of-fit diagnostics for regression models[END_REF], [START_REF] Staniswalis | Diagnostics for assessing regression models[END_REF]. However, these methods present several drawbacks. Some of them require some a priori choice of the regularity of µ, and the result of the test depends on this choice. Moreover, a data driven choice of the smoothing parameter involved in the non-parametric estimation, like the bandwidth in a kernel estimator, for example, may affect the level of the test. Another procedure proposed by [START_REF] Dette | Validation of linear regression models[END_REF] is based on the estimation of the empirical L 2 distance between µ and the null hypothesis, and does not depend on the choice of such a smoothing parameter. [START_REF] Horowitz | An adaptive, rate-optimal test of a parametric mean-regression model against a nonparametric alternative[END_REF] proposed a procedure that rejects the null hypothesis if the distance between the nonparametric kernel estimator and the kernel-smoothed parametric estimator of µ under the null hypothesis is large for some bandwidth within a grid. This test asymptotically achieves the desired level. It is rate-optimal among adaptive procedures over Hölder classes of alternatives and almost achieves the parametric rate of testing for directional alternatives. A similar approach was proposed by [START_REF] Guerre | Data-driven rate-optimal specification testing in regression models[END_REF].

An alternative approach is based on the cumulative sums of residuals over covariates or predicted values of the response variable. Due to their construction, such procedures do not depend on any smoothing parameter. They have shown themselves to be of an asymptotic level equal to the nominal level. They are consistent and optimal for testing the null hypothesis versus directional alternatives at the rate of n -1/2 ( see [START_REF] Diebolt | On testing the goodness-of-fit of nonlinear heteroscedastic regression models[END_REF], [START_REF] Stute | Nonparametric model checks for regression[END_REF], [START_REF] Su | A lack-of-fit test for the mean function in a generalized linear model[END_REF]).

A different approach proposed by [START_REF] Crainiceanu | Likelihood ratio tests in linear mixed models with one variance component[END_REF] for testing a polynomial function versus a parametric alternative consists in testing for the presence of a random effect, the alternative being characterized by splines with random coefficients. This approach requires the choice of the degree and knots of the splines. The test procedure is based on likelihood ratio test statistics and on their exact distributions under the null hypothesis. A similar procedure based on wavelets was proposed by [START_REF] Claeskens | Lack-of-fit tests in linear mixed models with application to wavelet tests[END_REF]. [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] proposed a multiple testing procedure for the expectation of the response variable Y in the case of a homoscedastic one-dimensional Gaussian regression with an unknown variance parameter. It is based on a large collection of Fisher tests, each of them testing that f belongs to V against a parametric alternative of possibly high dimension. The size of the test is exactly equal to the nominal level. The authors showed that the procedure is simultaneously consistent and rate-optimal over various classes of alternatives.

In the case of mixed-effects models, many papers have dealt with the problem of testing the random effects, but only a few papers have focused on the problem of testing the fixed effects. A review of methods for testing polynomial covariate effects in linear and generalized linear mixed models is given by [START_REF] Huang | Testing polynomial covariate effects in linear and generalized linear mixed models[END_REF].

In the case of generalized linear mixed-effects models, following the approach of [START_REF] Härdle | Testing parametric versus semiparametric modeling in generalized linear models[END_REF], [START_REF] Lombardía | Semiparametric inference in generalized mixed effects models[END_REF] proposed a procedure that combines the fully parametric likelihood approach for random effects models with a semiparametric regression profiled likelihood. The test statistic is based on the distance between estimates of the function µ under the null hypothesis and under the alternative. The asymptotic properties under the null hypothesis are established. A parametric bootstrap procedure is proposed for estimating the quantiles of the statistics distribution. [START_REF] Pan | Goodness-of-fit methods for generalized linear mixed models[END_REF] generalized the work of [START_REF] Su | A lack-of-fit test for the mean function in a generalized linear model[END_REF] and [START_REF] Stute | Nonparametric model checks for regression[END_REF], and proposed several procedures based on the cumulative sums of residuals. As in the case of homoscedastic regression, their procedures do not depend on the choice of a smoothing parameter. They show that their tests are asymptotically of the nominal level. In the particular case of linear mixed-effects models, they show the consistency of their procedure when the marginal mean is misspecified. [START_REF] Scheipl | Size and power of tests for a zero random effect variance or polynomial regression in additive and linear mixed models[END_REF] generalized the procedure proposed by [START_REF] Crainiceanu | Likelihood ratio tests in linear mixed models with one variance component[END_REF] for linear mixed-effects models. Following the work of [START_REF] Greven | Restricted likelihood ratio testing for zero variance components in linear mixed models[END_REF], they proposed to compute a pseudolikelihood by plugging an estimator of the random part of the mixed-effects model into the complete likelihood. To our knowledge, no theoretical property has been established in that context. Recently, [START_REF] Greven | On likelihood ratio testing for penalized splines[END_REF] established theoretical results for the score test statistic under the null hypothesis, the dimension of the spline basis being allowed to increase with the sample size. Nevertheless, the asymptotic distribution under the null hypothesis is not available in the case of mixed models.

Model and testing procedure

We consider the Gaussian regression model with block correlated errors. Let Y 1 , . . . , Y n be n Gaussian independent random vectors so that for each i = 1, . . . , n, E(Y i ) = f i and Var(Y i ) = Σ i (γ). The vectors (f i ) 1≤i≤n are unknown vectors of R J , and the matrices (Σ i ) 1≤i≤n are symmetric positive matrices that depend on q unknown parameters designated as γ. All the results presented in this paper remain true in the general case where the dimensions J i are not all equal to J as soon as max i J i is smaller than some constant.

The random vector Y satisfies:

Y = f + W 1/2 (γ)ε, (6.1)
where ε is a centered, standardized Gaussian vector of R nJ , and W = W (γ) is the block diagonal matrix with Σ i (γ) in block i. Our aim is to test the null hypothesis that f belongs to V , where V is a linear subspace of R nJ of dimension p. Many regression models in which it is of interest to test goodness-offit may satisfy Equation (6.1). Such models have been widely studied in the literature (e.g., [START_REF] Davidian | Nonlinear models for Repeated Measurement Data[END_REF], [START_REF] Demidenko | Mixed models. Wiley Series in Probability and Statistics[END_REF], [START_REF] Diggle | Analysis of longitudinal data[END_REF], [START_REF] Jones | Longitudinal data with serial correlation: a state-space approach[END_REF], [START_REF] Pinheiro | Mixed-effects models in S and S-Plus[END_REF], [START_REF] Seber | Nonlinear regression[END_REF], [START_REF] Vonesh | Linear and nonlinear models for the analysis of repeated measurements[END_REF]). Let us quote for example linear mixed effects models, hierarchical models, heteroscedastic regression models, regression models with autocorrelated errors, growth curves models.

We propose a goodness-of-fit test to test that the expectation of Y belongs to a specific linear subspace. We generalize the procedure developed by [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] to our framework as follows: we define a collection of alternative hypotheses whose cardinality may depend on n, and build on parametric models with low and high dimensions, diversified enough to cover a wide variety of possible alternatives. For each of these alternatives, we consider the likelihood ratio test statistic that would be obtained if the covariance matrix W were known. We then replace the unknown covariance matrix in this statistic by its maximum likelihood estimate. The distribution of each statistic under the null hypothesis is approximated by a χ 2 distribution. Using the Bonferroni adjustment, the null hypothesis is rejected as soon as it is rejected by one of the parametric tests.

The heuristic of our testing procedure can be described in the following way. Let us first assume that the matrix W is known, and consider the test of level α ∈]0, 1[ of the hypothesis H 0 : f ∈ V against the alternative H : f ∈ V + S where S is a linear subspace of R nJ such that the projection of S onto the space orthogonal to W -1/2 V in R nJ is non zero. The likelihood ratio test statistic (see [START_REF] With Discussion | Theoretical statistics[END_REF]) denoted by T S (W, Y ) is defined as:

T S (W, Y ) = Π S Π (W -1/2 V ) ⊥ W -1/2 Y 2 . The hypothesis H 0 is rejected if T S (W, Y ) > χ-1 D (α)
, where D is the dimension of Π (W -1/2 V ) ⊥ S and where χ-1 D (α) designates the (1 -α)-quantile of a χ 2 distribution with D degrees of freedom. Since the matrix W is unknown, we propose to plug the test statistic into the maximum likelihood estimator W of W under the alternative hypothesis H.

Following the idea proposed by [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] in the linear Gaussian regression model, the multiple testing procedure is constructed as follows. Let us consider a collection of linear subspaces of R nJ denoted as {S m , m ∈ M} where M is a set of indices that depends on n. We denote the hypothesis f ∈ V + S m by H m and the maximum likelihood estimator of W under H m by W m . We assume that for each m ∈ M, the dimension D m of Π ( W -1/2 m V ) ⊥ S m is greater than 1. Following the Bonferroni procedure, we choose some sequence {α m , m ∈ M} of positive numbers satisfying m∈M α m = α. The statistic for testing that f belongs to V against that it does not is then defined as:

T (α) = sup m∈M T m ( W m , Y ) -χ-1 Dm (α m ) , (6.2)
where T m stands for T Sm . The hypothesis H 0 is rejected when T (α) is positive.

Theoretical results

We evaluate the properties of the procedure when the sample size n tends to infinity, taking the fact that the cardinality of the collection as well as the dimensions of the alternatives are allowed to grow with the sample size into account. We show that the test is asymptotically of the desired level and that it is consistent over a large class of alternatives. In particular, the rates of testing are the same as those obtained in the homoscedastic case.

Asymptotic level and power

Consider the regression model as defined in Equation (6.1). Let V be a linear subspace of R nJ , and let {S m , m ∈ M} be a collection of linear subspaces as defined in Section 6.2. Let {α m , m ∈ M} be a sequence of positive numbers such that m∈M α m = α.

Assuming some regularity assumptions on the model and on the collection of linear subspaces, we have the following theorem: Theorem 6.1. Suppose that D 4 m /n tends to 0 and m (D m /n) 1/3 tends to 0. Consider the statistic T (α) as defined in Equation (6.2), then:

lim n→+∞ P H0 (T (α) > 0) ≤ α,
where P H0 denotes the probability when f ∈ V .

Following are some comments about the assumptions that make it possible to prove this theorem. Assuming that D 4 m /n tends to 0, it is possible to show the existence of the parameter estimators defined under H m . Moreover, if θ V denotes the p coefficients of Π V f in V , it can be shown that the estimator of the parameters θ V is n/D m -consistent. Because the modeling of the covariance matrix W does not depend on S m , and as a result of the block diagonal structure of W , the estimator of γ is √ n-consistent. Finally the statistic T m ( W m , Y) is asymptotically equal to a χ 2 variable with D m degrees of freedom up to a remainder term that is of order D m / √ n. This result follows from the consistency of the parameter estimators and from the regularity of the density of a χ 2 variable. Assuming that m (D m /n) 1/3 tends to 0, it is possible to control the discrepancies between T m ( W m , Y) and the χ 2 over all alternatives H m .

Let us now consider the power of the test. Assuming some regularity assumptions on the model and on the collection of linear subspaces , we have the following result: Theorem 6.2. Suppose that for all m ∈ M, D 4 m /n tends to 0. Let E n be the set of f ∈ R nJ for which there exists m ∈ M such that:

1 n f -Π V +Sm f 2 = O 1 √ n , (6.3) 1 n Π (W -1/2 V ) ⊥ W -1/2 f 2 ≥ 1 n Π S ⊥ m Π (W -1/2 V ) ⊥ W -1/2 f 2 + v 2 n,m , with v 2 n,m = κ n D m log log(n) α m + log log(n) α m (1 + o(1)) ,
for some positive constant κ. Then:

lim n→+∞ sup f ∈En P f (T (α) ≤ 0) = 0,
where P f denotes the probability under Model (6.1).

Following are some comments about the theorem. For each m ∈ M, under the assumption that D 4 m /n tends to 0, the estimators of the parameters under H m converge in probability under P f . Their rates of convergence are the same under P f as under P H0 .

The condition f -

Π V +Sm f 2 = O ( √ n) results from remainder terms in the asymptotic expansion of T m ( W m , Y) of the following form: Π Sm Π (W -1/2 V ) ⊥ ( W -1/2 m -W -1/2 )f 2 .
It can be shown that these terms are of the order f -Π V +Sm f / √ n. If this quantity does not vanish, then we are not able to give an expansion for T m ( W m , Y) under P f .

Finally, it should be emphasized that the power of our procedure is comparable to the power of the exact procedure proposed by [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] in the case of Gaussian linear regression variance with independent and homoscedastic errors.

Detection of local alternatives

The aim of this section is to establish the power of the test for vectors f such that f -Π V f is of order √ log log n.

We consider a linear space S such that the dimension D of Π (W -1/2 V ) ⊥ S is greater than 1, and we denote the maximum likelihood estimator of W under the assumption f ∈ V + S by W .

Let us denote the canonical basis of R nJ by {e i,a , a = 1, . . . , J, i = 1, . . . , n}. The collection {S m , m ∈ M} is composed here of all vectors of the canonical basis such that the dimension D m of Π ( W -1/2 V ) ⊥ S m is equal to 1. The collection {α m , m ∈ M} is defined as before.

Finally, we consider the test statistic:

U (α) = sup m∈M T m ( W , Y ) -χ-1 1 (α m ) . (6.4)
Assuming some regularity assumptions on the model and on the collection of linear subspaces, we have the following result. Theorem 6.3. Suppose that D 4 /n tends to 0 and that for all m ∈ M, α m ≥ exp(-n/(log n) 3 ). If we consider the statistic U (α) as defined in Equation (6.4), then:

lim n→+∞ P H0 (U (α) > 0) ≤ α,
where P H0 denotes the probability when f ∈ V .

For g ∈ R nJ such that g 2 /n is bounded, let:

A n (g) = f ∈ R nJ such that f = Π V f + log log n n g .
If f ∈ A n (g), and if there exists m 0 such that:

Π Sm 0 Π (W -1/2 V ) ⊥ W -1/2 g 2 ≥ Cn log 1 α m0
, (6.5)

then:

lim n→+∞ P f (U (α) ≤ 0) = 0 ,
where P f denotes the probability under Model (6.1).

The Bootstrap procedure

We also propose a bootstrap procedure for estimating the quantiles of the distribution of each test statistic and prove that this bootstrap procedure is asymptotically of the desired level.

In the preceding sections, we proposed a test procedure based on the approximation of the quantiles of T m ( W m , Y) by those of a χ 2 variable with D m degrees of freedom. An alternative to this procedure is the bootstrap, where the quantiles of T m ( W m , Y) are approximated by the quantiles of their bootstrap distribution under H 0 .

Let ε be a centered, standardized Gaussian vector of R nJ independent of ε, and let:

Y = Π V Y + W 1/2 0 ε , (6.6)
where W 0 is the maximum likelihood estimator of W under H 0 .

Let

T m = T m ( W m , Y )
, where W m is the bootstrap version of W m , and let q m (α) be the (1-α)-quantile of the distribution of T m conditional on Y . We then consider the test statistic defined as:

T (α) = sup m∈M T m ( W m , Y ) -q m (α m ) .
(6.7)

The hypothesis H 0 is rejected when T (α) is positive.

The proof that the asymptotic level of the bootstrap procedure is α is similar to the proof of Theorem 6.1. However, it needs to uniformly control the difference between the distribution function of T m and its bootstrap version. This leads to an additional assumption about the collection M, precisely, about the cardinality of the subset M 1 defined as the set of m ∈ M such that D m = 1. Theorem 6.4. Assume that the assumptions of Theorem 6.1 are fulfilled and that |M 1 |/n 1/5 tends to 0. Let the statistic T (α) be defined as in Equation (6.7), then:

lim n→+∞ P H0 (T (α) > 0) ≤ α.

Numerical studies

We conducted a simulation study to to assess the performances of our procedure when n is fixed, and to compare it with the omnibus test proposed by [START_REF] Pan | Goodness-of-fit methods for generalized linear mixed models[END_REF] based on cumulative residuals. We consider the case of a functional regression model with random effects where we observe (Y ij , X ij ) for i = 1, . . . n and j = 1, . . . J under the following model:

Y ij = µ(X ij ) + b i + ε ij ,
where the variables b i and ε ij are independent centered Gaussian variables with the variance denoted as σ 2 b and σ 2 , respectively. Our aim is to test that µ is constant. In that case, V is the linear space spanned by the column vector with all components equal to 1.

We consider three types of alternatives, the first with smooth variations, the second with heavy variations and the third with oscillations. We study the level and the power of the test, as well as the effect of the choice of the collection of linear subspaces.

The simulation study shows that the bootstrap procedure is particularly efficient in the case of a small sample size. Our tests give results on the same order as those of Pan and Lin in the case of regular deviations from the null hypothesis and outperform their method in the context of oscillating deviations from the null hypothesis.

We complete our experiments by considering a real dataset of forest coverage in Galicia already treated in [START_REF] Lombardía | Semiparametric inference in generalized mixed effects models[END_REF]. This application highlights, in particular, the advantage of our procedure of being independent of the choice of the smoothing parameter.

Perspectives

It would be interesting to investigate further works in this area. First, one could consider a non linear hypothesis instead of a linear one. The test procedure could be extended possibly using an approximation function basis. Second, one could investigate the case where the size J of each Gaussian vector goes to infinity. This would lead to test functional hypothesis.

Part III

Modeling and Inference in Survival Analysis

In the third part, I included my contributions dealing with modeling and estimation in survival analysis. Related studies include [A6,A7,P6,R1].

In collaboration with Charles El-Nouty, we considered the frailty models introduced by [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF] which are an extension of the Cox model that takes the heterogeneity that exists in survival data into account by introducing latent variables. We were interested in these models since we have been in contact with researchers at the Unité de Recherche en Epidémiologie Nutritionnelle (UREN) at the University of Paris 13, who often need to analyze this type of dataset. Very rich from a modeling point of view, these models are very complex from a mathematical point of view, and the estimation task is often very difficult. In numerous cases, the existing algorithms do not provide satisfactory solutions. We proposed the application of the SAEM-MCMC algorithm to frailty models to evaluate the MLE [A6]. We proved that under general regularity assumptions fulfilled by classical frailty models, the SAEM-MCMC algorithm is almost surely convergent toward a local maximum of the observed likelihood. We compared the performances of this algorithm with others that exist in the literature on simulated data and on a real set of bladder cancer data. The numerical results showed a net advantage when using the SAEM-MCMC algorithm, both in terms of the accuracy of the limit as well as the computation time.

In collaboration with Luc Duchateau and Klaartje Goethals (Faculty of Veterinary Medicine, Ghent University), we have analyzed a mastitis epidemic dataset using frailty models with a frailty vector of size four with different covariance structures [R1]. Assessing the correlation structure in cow udder quarter infection times allows us to analyze the propagation risk of the disease as a function of the position of the infection. The parameter estimation that could not be performed in a reasonable time in such models until now, is done using the SAEM-MCMC algorithm. We compared four nested models using likelihood ratio tests. Using simulation studies, we justified a posteriori the use of likelihood ratio tests for finite sample size.

Some survival models

Context and notation

Survival analysis consists mainly in the study of event times of interest for different individuals of a given population. For example, this could be failure times such as the death of some patients from a given disease in epidemiology, the breakdown of some systems in reliability, or the appearance of some particular leaves of interest for a plant in agronomy. The analysis of event times attempts to answer many questions such as: which proportion of a population will survive past a certain time? How will the surviving individuals die or fail? Are there particular characteristics that increase or decrease the chance of survival? Are there multiple causes of death or failure?

Let us introduce some general notations useful in the context of survival analysis. We consider a population of N individuals. We denote by T i the random variable corresponding to the event time of interest, also referred to as lifetime, of the i-th individual of the population for 1 ≤ i ≤ N .

A crucial quantity of interest is the survival function of the i-th individual defined for t > 0 by S i (t) = P (T i ≥ t). The most popular estimator of the survival function based on the lifetime data is the Kaplan-Meier estimator (see [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF]) that has been studied by many authors.

A second quantity of interest is the hazard rate λ i of the i-th individual defined for t > 0 by:

λ i (t) = lim dt→0 + P (t ≤ T i < t + dt|T i ≥ t) dt ,
which is the probability for the i-th individual that the event of interest arises at time t given that it has not arisen before t. This quantity can be interpreted in epidemiology, for example, as the risk of dying from a given disease, a constant hazard rate corresponding to a chronic disease.

The survival function can be related to the hazard rate as follows. Let us assume that the random variable T i has a probability density function denoted by f i . Then, for t > 0, we have:

λ i (t) = lim dt→0 + P (t ≤ T i < t + dt|T i ≥ t) dt = 1 S i (t) lim dt→0 + S i (t) -S i (t + dt) dt = f i (t) S i (t) = -[log S i (t)] ,
thus

S i (t) = exp - t 0 λ i (s)ds .
Note that a constant hazard rate equal to λ will correspond to an exponential distribution of parameter λ for the survival time T i .

As mentioned before, it could be interesting to link some particular characteristics of the individuals of the population to their behavior when faced with the event of interest. For example if we consider several treatments for a population of patients, then covariates such as the type of treatment, age and gender of the patient can be taken into account in the model. Some survival models make it possible to include covariates in the analysis of lifetimes. The most popular one is the Cox model (see [START_REF] Cox | Regression models and life-tables[END_REF]).

The Cox Proportional Hazards Model

The Cox Model, also known as the Proportional Hazards Model, states that covariates act multiplicatively on the hazard rate. Thus, the hazard rate is explained as follows:

λ i (t) = λ 0 (t) exp(x i β), (7.1)
where X i are the covariates of individual i, λ 0 is an unknown baseline function, and β is a parameter effect vector. Note that two individuals sharing the same covariates have the same hazard rate. However, the lifetimes of two individuals sharing the same covariates differ due to the hazard of the lifetime distribution.

The main assumption of the Cox model is the so-called risk proportional assumption.

Indeed, as a consequence of Equation (7.1), we have: Thus, if we consider two groups A and B of individuals with the same covariates within each group, the plot of the function t → log(-log(S A (t)) should be obtained as the vertical translation of the one of t → log(-log(S B (t)). This is illustrated on simulated data in Figure 12. However, the datasets do not often satisfy this assumption. For example, Figure 13 presents the plot of the survival function of real data that do not satisfy this assumption. A classical underlying assumption of the Cox model is that the observations are independent, at least conditionally on covariates. However, this assumption is often not fulfilled by the data because of the lack of homogeneity for the population of interest. For example it can happen that failure times are clustered into groups such as families or geographical areas. Illustrations on medical datasets are given in [START_REF] Aalen | Analysing incidence of testis cancer by means of a frailty model[END_REF] where the authors studied the incidence of age on testiculars cancer, or in [START_REF] Gray | Spline-based tests in survival analysis[END_REF][START_REF] Gray | Tests for variation over groups in survival data[END_REF] where the author showed the effects of participating institutions in a multicenter lung cancer trial. The same frameworks are developed in economy (see [START_REF] Horowitz | Semiparametric and nonparametric methods in econometrics[END_REF]).

log S i (t) = -
To face this lack of fit of the Cox model, one solution is to take the heterogeneity present in the data into account by introducing some random effects in the modeling of the hazard rate. This is the aim of frailty models introduced by [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF].

The Frailty Model

Let us introduce some notations used in the context of frailty models. We consider a population of individuals clustered into N groups. For 1 ≤ i ≤ N , we denote by n i the size of the i-th group. For 1 ≤ i ≤ N and 1 ≤ j ≤ n i , the event time for the subject j of the group i is modeled by a random variable denoted by T ij . We also consider here the censoring time for the subject j of the group i which is modeled by a random variable denoted by C ij . However, the event time (T i ) and the censored time (C i ) are generally not observed. Let us also define the random variables Y ij = min(T ij , C ij ) and ∆ ij = 1 {Tij ≤Cij } , where 1 A denotes the indicator function of any set A. We observe the couples (Y ij , ∆ ij ) 1≤i≤N,1≤j≤ni . Let b i be the common frailty random vector for the i-th group for 1 ≤ i ≤ N . The initial choice of [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF] of the Gamma distribution for the frailty was motivated by its mathematical convenience. It was extended by [START_REF] Clayton | Multivariate generalizations of the proportional hazards model[END_REF] to other frailty distributions. We refer to [START_REF] Hougaard | Analysis of multivariate survival data[END_REF] and to [START_REF] Duchateau | The Frailty Model[END_REF] for further information on the choice of the frailty law.

Denote by β an unknown parameter effects vector and by λ 0 the unspecified baseline hazard function. Let x ij and z ij be design vectors of covariates associated with failure. Consider λ ij (t|b i ) the conditional hazard function for the j-th individual of the i-th group at time t. The traditional frailty model is defined by:

λ ij (t|b i ) = λ 0 (t) exp(x t ij β + z t ij b i ), (7.2)
for 1 ≤ i ≤ N, 1 ≤ j ≤ n i and t ≥ 0, where t denotes the transposition operator.

The classical assumptions on the frailty model are given below:

• F1 : The censoring times (C ij ) 1≤i≤N,1≤j≤ni are independent of the event times (T ij ) 1≤i≤N,1≤j≤ni and of the frailties (b i ) 1≤i≤N . • F2 : Conditional to the frailties (b i ) 1≤i≤N , the event times (T ij ) 1≤i≤N,1≤j≤ni are independent.

• F3 : The frailties (b i ) 1≤i≤N are independent and identically distributed with probability density function f η on R q , where η is an unknown vector.

We denote the vector of all unknown parameters by θ = (α, β, η).

Recall that the frailties are unobserved random variables. The frailty model (7.2) can therefore be seen as a linear mixed effects model for the logarithm of the hazard function. In the terminology of mixed effects models, the frailties (b i ) 1≤i≤N correspond to random effects and the vector β to fixed effects. Thus, the frailty model can also be viewed as a proportional hazards model with random effects, also referred to as proportional hazards mixed model (see [START_REF] Vaida | Proportional hazards model with random effects[END_REF]).

Estimation in Frailty Model

We choose here a frequentist approach and consider the MLE for the parameter θ, namely the value θN of θ which maximizes the marginal likelihood denoted by L obs N (see [START_REF] Duchateau | The Frailty Model[END_REF]). This quantity is obtained by integrating the complete likelihood L N over the unobserved frailties, given by:

L N (y, δ, b; θ) = N i=1 f η (b i ) × N i=1 ni j=1 λ ij (y ij |b i ) δij exp - yij 0 λ ij (u|b i )du , (7.3)
where we denote the vectors corresponding to the realizations (y ij ), (δ ij ) and (b i ) by y, δ and b, respectively. Asymptotic theoretical properties for the Maximum Likelihood Estimator (MLE) were established by [START_REF] Murphy | Consistency in a proportional hazards model incorporating a random effect[END_REF], [START_REF] Murphy | Asymptotic theory for the frailty model[END_REF] and [START_REF] Parner | Asymptotic theory for the correlated gamma-frailty model[END_REF]. The authors proved under general conditions that θN exists and converges toward θ as N goes to infinity with probability one. Nevertheless it is often not possible in many practical cases to compute the MLE directly. To overcome this crucial difficulty, two main approaches were developed.

The first one consists in applying Cox's idea to the frailty models in order to obtain an approximated marginal likelihood. Examples where this approach is used can be found in [START_REF] Mcgilchrist | Regression with frailty in survival analysis[END_REF] (penalized likelihood), in [START_REF] Nielsen | A counting process approach to maximum likelihood estimation in frailty models[END_REF] (partial likelihood), in [START_REF] Therneau | Modeling Survival Data : Extending the Cox Model[END_REF] (penalized partial likelihood) and in [START_REF] Rondeau | Maximum penalized likelihood estimation in a gamma-frailty model[END_REF] (penalized full likelihood). This approach has been extended to a Bayesian model introduced by [START_REF] Ducrocq | A bayesian analysis of mixed survival models[END_REF] and developed by [START_REF] Legrand | A Bayesian approach to jointly estimate centre and treatment by centre heterogeneity in a proportional hazards model[END_REF] and [START_REF] Legrand | Validation of prognostic indices using the frailty model[END_REF].

Consider now the second approach. It consists in a numerical approximation of the MLE instead of a direct computation. Since the frailties are not observed, the underlying model belongs to the family of models with hidden variables. Thus, a powerful tool to solve the MLE problem is the EM algorithm proposed by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. However, in many frailty models, the EM algorithm cannot be directly applied, more specifically, the expectation step (E-step). Thus, several authors suggest approximations of this algorithm. We mention two of them here that are commonly used. [START_REF] Ripatti | Maximum likelihood inference for multivariate frailty models using an automated Monte Carlo EM algorithm[END_REF] applied the Monte Carlo EM (MCEM) algorithm to frailty models. It was introduced by [START_REF] Wei | A Monte Carlo implementation of the EM algorithm and the Poor's Man's data augmentation algorithms[END_REF], whereas its convergence in a general setting was established by [START_REF] Fort | Convergence of the Monte Carlo Expectation Maximization for curved exponential families[END_REF]. However, it requires intensive computation time. Cortiñas [START_REF] Cortiñas Abrahantes | A version of the EM algorithm for proportional hazard model with random effects[END_REF] then applied the deterministic EM-Laplace to frailty models. Unfortunately this method requires that the number of observations sharing the same frailty tends to infinity. In practice, this assumption is often not fulfilled, for example, in twin studies. Moreover, this method can induce some bias (see Cortiñas Abrahantes and Burzykowski (2005, pp. 853-859)). A well-known example of the lack of convergence is given in [START_REF] Allassonnière | Towards a coherent statistical framework for dense deformable template estimation[END_REF] (pp. 12-13). As a result, it must be used very carefully. As far as we know, there is no theoretical proof of the convergence of the other existing approximations.

Estimation in Frailty Model using the SAEM-MCMC Algorithm

Algorithmic method

In collaboration with Charles El-Nouty,we proposed to use the Stochastic Approximation Expectation Maximization with the Monte Carlo Markov Chain (SAEM-MCMC) algorithm introduced in (A1) for maximum likelihood estimation (MLE) in frailty models (see [A6]). We briefly recall here the characteristics of this method detailed in Section 2. The usual expectation step of the EM algorithm is divided into two new steps: the first one consists in simulating one realization of the non observed frailties, whereas the second one computes a stochastic approximation of the complete log-likelihood by using this simulated value of the frailties. The maximization step follows the same lines as those of the EM algorithm.

Moreover, this algorithm of convenient use has theoretical properties and requires less simulations than the MCEM algorithm. The proposed algorithm also makes it possible to deal with the multivariate frailty models without assumption on the frailty covariance structure. Another advantage of the SAEM-MCMC algorithm is that it generates the estimation of the observed Fisher information matrix at the same time.

Convergence property of the algorithm

I detail here the assumption required to ensure the convergence of the proposed estimation algorithm, as well as the related convergence result.

Let us introduce two additional assumptions on the frailty model denoted (F4-F5). The first concerns the baseline function λ 0 , whereas the second one deals with the frailty distribution f η . Extra assumptions denoted by (H1-H2) and by (SAEM1-SAEM3) in [A6] are made on the regularity of the model and on the stochastic approximation procedure, respectively (see [A6] for more details). Our convergence result is given in the following theorem. Note that under additional regularity assumptions, the almost sure convergence of the sequence (θ k ) k toward a local maximum of log L obs N is obtained (see [A1]). The usual choices made for the parametric baseline λ 0 (respectively for the frailty distribution) satisfy the assumption (F4) (respectively (F5)). We refer to [START_REF] Duchateau | The Frailty Model[END_REF] for further information about this field. For example, when the frailty follows a Gaussian, a Gamma or a Weibull distribution, the SAEM-MCMC algorithm converges almost surely toward the MLE under weak additional regularity conditions.

Experiments on a bladder cancer dataset and simulation studies

Experiments a the bladder cancer dataset

Let us briefly describe the bladder cancer dataset (EORTC trials 30781,30782,30791,30831,30832,30845 and 30863 Genito-Urinary Tract Cancer Group) that we deal with. It is composed of 2596 eligible patients recruited by 63 medical centers. One hundred patients having missing values and 24 centers have less than 20 patients. In our analysis, we keep 2286 patients and 39 medical centers with more than 20 patients as in [START_REF] Abrahantes | Comparison of different estimation procedures for proportional hazards model with random effects[END_REF] and in [START_REF] Legrand | A Bayesian approach to jointly estimate centre and treatment by centre heterogeneity in a proportional hazards model[END_REF].

The study of the dataset shows that it is a setting with approximately 51% of censuring and that approximately 80% of the individuals follow an intravesical treatment (see [START_REF] Sylvester | Predicting recurrence and progression in individual patients with stage ta t1 bladder cancer using eortc risk tables: A combined analysis of 2596 patients from seven eortc trials[END_REF]). To obtain a satisfactory Laplace approximation, we have to choose a very sharp covering (or tolerance) of order 10 -6 to find a global minimum and to assume regularity conditions on the conditional distribution of the frailty variables. Otherwise, many trajectories do not converge. Moreover, the nature of the trajectory is deeply affected by small variations of the initial values. This phenomenon is not really surprising and is also highlighted in [A3] where the EM-Laplace algorithm is termed as FAM-EM algorithm. Consider now the computing time problem. As soon as the EM-Laplace converges, it requires the same computation time as the SAEM-MCMC algorithm. However, when the EM-Laplace does not converge, the computation time can increase considerably. Finally, the use of the MCEM algorithm requires a computation time ten times greater than the SAEM-MCMC, taking the number of Monte Carlo simulations equal to 5 during the first 50 iterations and equal to 1000 thereafter. This is the main difficulty to apply the MCEM algorithm and therefore justifies the use of the SAEM-MCMC algorithm.

Let us now comment on the numerical results. As expected, the two stochastic algorithms give same order estimates of parameters and of the marginal log-likelihood, whereas those obtained by using the deterministic algorithm are very different. The EM-Laplace algorithm leads to an estimation of the marginal log-likelihood by -9490.47, which is strictly lower than the estimations of the marginal log-likelihood obtained with the MCEM algorithm and with the SAEM-MCMC algorithm, -9366.59 and -9365.41, respectively. Note that these two last values are local maxima of the marginal loglikelihood. The estimates of the parameters (10.563 × 10 -4 , -0.0167, 0.4074, 0.0032) given by the EM-Laplace algorithm do not maximize the marginal loglikelihood and are consequently meaningless. The estimated standard errors given by the two stochastic algorithms have the same order for λ 0 and β, whereas there is a small gain for σ 2 0 and σ 2 1 when using the MCEM algorithm. Some studies for the same real dataset suggest that for any i = 1, . . . , 39, the center effect b 0i could be correlated with the treatment per center interaction b 1i . For this purpose, we consider a second model, referred to as Model II and defined as follows:

λ ij (t|b i ) = λ 0 (t) exp b 0i + x t ij (β + b 1i ) ,
for 1 ≤ i ≤ N, 1 ≤ j ≤ n i and t ≥ 0, where the N random vectors (b t i ) 1≤i≤N are independent and identically distributed multivariate Gaussian:

b i = b 0i b 1i ∼ N 2 0 0 , σ 2 0 σ 01 σ 01 σ 2 1
We insist on the fact that we introduce a dependence structure on the two frailty terms. To study Model II, we limit ourselves to the SAEM-MCMC algorithm. In Table 5 we recall the results obtained for Model I in Table 4 and give the new numerical results for Model II. Note that the parameter σ 01 only exists in Model II.

We would like to make some comments on these results at this point. Models I and II have the same order estimates for λ 0 and β, whereas the estimates of σ 2 0 and σ 2 1 are lower when considering Model II. This is not really surprising because there is a covariance term σ 01 . Concerning standard error terms, there is a big gain when considering Model II, particularly for the parameters σ 2 0 and σ 2 1 . Based on the numerical results obtained for Model II, it is natural to determine if the covariance σ 01 is significantly different from zero. Since Model I and Model II are nested, we do a likelihood ratio test where the null hypothesis is "σ 01 = 0". The numerical value of the test statistic is -2 -9365.41 + 9357.51 = 15.80. Let Z be a χ 2 (1) random variable. Since 15.80 > 6.63 where P (Z > 6.63) = 0.01, we reject the null hypothesis at level 1%. Finally, we recommend the use of Model II to analyze this dataset.

Simulation studies

We conduct a simulation study to highlight the performance of our method regarding the other ones that exist in the literature in the same setting as the analysis of the bladder cancer dataset carried out above.

We consider two different settings of censoring namely a moderate one and a heavy one. The design vectors x ij are equal to 1 for 70% of the individuals and 0 for 30% of the individuals in the moderate censoring setting with 40% of censoring, and equal to 1 for 50% of the individuals and 0 for 50% of the individuals in the heavy censoring setting with 60% of censoring. Although the parameter λ 0 could be estimated, we set λ 0 = 0.077 in order to keep the same choice as the one made in [START_REF] Abrahantes | Comparison of different estimation procedures for proportional hazards model with random effects[END_REF]. We generate R = 250 datasets for each setting. Let us denote by θ (r) = (0.077, β (r) , σ 2(r) 0 , σ

2(r) 1

) the estimates obtained on the r-th simulated dataset for 1 ≤ r ≤ R. Hence, the estimate θ is the empirical mean of the (θ (r) ) 1≤r≤R . The standard error of θ is computed by two methods. The first one is the empirical standard error σ( θ) of this estimate whereas the second one is the so-called mean model-based standard error that corresponds to the square root of the empirical mean of the variance estimate given by the inversion of the Fisher Information Matrix. Numerical results are given in Tables 6 and7. For each algorithm and each parameter, we present three numbers: the first one corresponds to the estimation of the parameter, whereas the two last numbers correspond to the estimation of the standard error, computed by the two methods described above.

In both censoring settings the SAEM-MCMC algorithm provides a better numerical approximation of the true parameter value than the EM-Laplace and the MCEM algorithms in terms of bias and empirical standard error estimates. However, although the EM-Laplace algorithm often leads to smaller mean model-based standard errors estimates, it also leads to bigger differences between the empirical and the mean model-based estimates of the standard errors. Note also that we obtain good approximations using the MCEM algorithm. However, to reach the same accuracy as in the results obtained by the SAEM-MCMC algorithm, we have to consider large sample sizes in the Monte Carlo approximation that lead to longer computation times.

Assessing the correlation structure in cow udder quarter infection times through extensions of the correlated frailty model

In collaboration with Charles El-Nouty (LAGA, University Paris 13 Sorbonne Paris Cité), Luc Duchateau and Klaartje Goethals (Faculty of Veterinary Medicine, Ghent University), we propose to assess the correlation structure in cow udder quarter infection times through several correlated frailty models [R1].

The main contribution of this work is to consider frailty models with an unknown correlation structures on the frailty vector of size 4 which describe precisely the random effects of each of the 4 udder quarters.

The usually correlated frailty models were only handling frailty vectors of size 2, because the parameter estimation could not be carry out efficiently. By using the convergent stochastic estimation algorithm presented in Section 8, the estimation task can now be achieve precisely in a reasonable time. We consider different possible correlation structures between all the frailty terms and analyze the performance of each corresponding model. We also compare these nested correlated frailty models by using the likelihood ratio test. We evaluate the performance of the estimation algorithm and of the finite sample size property of the likelihood ratio test on simulated data in this setting.

The existing correlated frailty model methodology is extended to clusters of size four and two different correlation structures. Wienke (2011) gives an excellent overview of the correlated frailty model. Most of the research on the correlated frailty model has been done for bivariate survival data, i.e., clusters of size two. Parameter estimates are in most cases obtained by rewriting the correlated frailty model in copula form, next estimate the marginal survival functions, and finally estimating the association parameter(s) after imputation of the marginal survival functions in the copula form. The cluster size of four allows us to investigate many more different correlation structures, and compare them to choose the most likely correlation structure generating the data , using the log likelihood ratio test. Furthermore, the used methodology, the SAEM-MCMC algorithm, is much more linked to the frailty model framework, and does not require the transformation to the copula format.

Mastitis dataset: times to infection in four clustered udder quarters

Mastitis or udder infection is economically the most important dairy cow disease in the western world (Seegers et al., 2003). Infections occur at udder quarter level, as the four udder quarters are fully separated, and are infected individually. From a point of view of controlling the disease in a herd, it is essential to know whether udder quarters have a higher risk of infection when one of the udder quarters is infected. With one of the udder quarters infected, it would be helpful to know as well whether each of the three other udder quarters have the same possibly increased risk, or whether the risk of infection differs according to the location of the udder quarter relative to the infected udder quarter.

In order to study such relationships between the four udder quarters, the individual udder quarters were followed up for one lactation period for infection, which generated clustered censored infection time data (Laevens et al., 1997). Obviously, the four udder quarters are clustered within a cow. Such multivariate data can be modeled in different ways to cope with the clustering structure [START_REF] Duchateau | The Frailty Model[END_REF], but in our modeling approach it is important to capture and quantify the association between the infection times of the udder quarters within a cow. The copula model and the frailty model are two possible approaches.

The mastitis dataset consists in 1196 cows. The udder quarters of each cow are followed up individually for time to infection in a lactation period. The udder quarter in which no infection occured are right censored at the time of the end of the lactation period. Two different covariates are introduced in the frailty model. Parity is considered at two levels, primiparous and multiparous, and is a cow characteristic. Location is either front or rear, and changes obviously from one quarter to another within a cow.

Several correlation structures for the frailty model

We consider general frailty models as defined in Section 7.3 whose conditional hazard rate satisfayes Equation (7.2) to analyze this dataset.

Our purpose is to focus on the effect of the relative location of the infected udder quarter on the propagation of the infection. Thus we are interested in frailty models having a frailty vector of size 4 having a specific correlation structure on its components.

Let us denote the 4-dimensional frailty random vector by b

i = (b i1 , b i2 , b i3 , b i4
) for the i-th cow for 1 ≤ i ≤ N and the vector of all frailty terms by b = (b i ) 1≤i≤N . We make the classical independence assumptions for the frailty models: the censoring times (C ij ) 1≤i≤N,1≤j≤4 are independent of the event times (T ij ) 1≤i≤N,1≤j≤4 and of the frailties (b i ) 1≤i≤N ; conditionally to the frailties (b i ) 1≤i≤N , the event times (T ij ) 1≤i≤N,1≤j≤4 are independent.

Denote by λ 0 (t) the unspecified baseline hazard function at time t. In our setting the constant hazard rate is not realistic for describing the time to udder quarter infection, as it is mentioned in [START_REF] Goethals | Modeling intervalcensored, clustered cow udder quarter infection times through the shared gamma frailty model[END_REF]. Therefore we assume that the baseline function has a Weibull parametric expression given by λ 0 (t) = λ 0 γt γ-1 for λ 0 > 0 and γ > 0.

Denote by x ij the 2-dimensional design vector of covariates associated with failure. This corresponds to the primiparous or multiparous status of the cow and the location of the udder quarter, namely front or rear. Let us introduce β = (β 1 , β 2 ) an unknown 2-dimensional vector corresponding to the two covariate effects. Note that the value of one covariate (location) changes within a cow whereas the other one (parity) changes between cows, i.e. the four udder quarters of a cow have the same parity.

Let us denote the conditional hazard function by λ ij (t|b i ) for the j-th individual of the i-th group at time t and the transposition operator by t . The frailty model M 1 , for 1 ≤ i ≤ N, 1 ≤ j ≤ 4 and t ≥ 0 is given by: λ ij (t|b i ) = λ 0 γt γ-1 exp(x t ij β + b ij ), (9.1) where the frailties (b i ) 1≤i≤N are independent and identically multivariate Gaussian distributed N (0, Σ) with positive definite covariance matrix Σ given by:

Σ = σ 2      1 ρ 1 ρ 2 ρ 3 ρ 1 1 ρ 3 ρ 2 ρ 2 ρ 3 1 ρ 1 ρ 3 ρ 2 ρ 1 1     
with σ 2 ≥ 0 and (ρ 1 , ρ 2 , ρ 3 ) ∈ [-1, 1] 3 . This leads to the following explicit restrictions on the parameters:

     |ρ 2 -ρ 3 | < 1 -ρ 1 |ρ 2 + ρ 3 | < 1 + ρ 1 0 < σ 2 .
We insist on the fact that the rank of Σ equals 4. This ensures that the Gaussian law of the frailty random vector is non degenerate, i.e., its support is equal to the whole space R 4 .

Θ = R + * × R + * × R 2 × R + * × [-1, 1] 3 .
Figure 15 shows the correlation structure between the random effects of the udder quarters corresponding to the covariance matrix Σ defined in the model M 1 . The correlation between the random effects of the front left and right quarters and the rear left and right quarters is denoted by ρ 1 . The correlation between the random effects of the front left and rear right quarters and the front right and rear left quarters is denoted by ρ 2 . The correlation between the random effects of the left front and rear quarters and the right front and rear quarters of the front left udder quarter is denoted by ρ 3 .

We focus now on three particular models derived from the correlated frailty model M 1 , by considering more and more specific correlation structures imposed on the frailty random vector. Each of them corresponds to a practical situation for the interactions between udder quarters.

We define the sub-model M 2 of model M 1 by putting the correlation ρ 3 equal to the correlation ρ 2 . The matrix Σ being positive definite, it follows:

     |ρ 1 | < 1 |ρ 2 | < 1+ρ1 2 0 < σ 2
.

We define the sub-model M 3 of model M 1 by putting the correlations ρ 2 and ρ 3 all equal to the correlation ρ 1 . Thus in model M 3 the correlation between each pair of random effects is equal to ρ 1 . Note that the model M 3 is also a sub-model of model M 2 .

Table 8

Parameter estimates and mean model based standard errors in parenthesis for the mastitis dataset for each of the four models M 1 , M 2 , M 3 and M 4 with λ 0 and γ the parameters of the Weibull baseline hazard, β 1 the effect of parity, β 2 the effect of location and σ 2 , ρ 1 , ρ 2 and ρ 3 the parameters of the random effect distribution. The last line gives the log-likelihood. Let us make some comments on the numerical results obtained from the real data set of mastitis. The parameter estimates of the baseline hazard function and of the covariates effects are similar, except for those obtained for model M 4 which has a very elementary frailty correlation structure. The residual variance parameter σ 2 is higher in the models M 1 , M 2 and M 3 , compare to model M 4 . This can be explain by the fact that this term grows to allow the model to fit the data under more constrain correlation structure. The correlation parameters ρ 1 and ρ 2 are stable in these three models and very close to 1. The estimated log-likelihoods are also similar except for the model M 4 which has a much lower one. and lead us to conclude via a likelihood ratio test that the model M 4 is very too simple for fitting this data set. The three other models give similar estimation results and log-likelihood values. We next compare model M 3 with model M 2 . The likelihood ratio test statistic equals 4.23, leading to the p-value P (X > 4.23) = 0.04, where X is a chi-squared random variable with one degree of freedom. Therefore, model M 3 can be rejected in favour of model M 2 at the 5% significance level. Nevertheless let us keep in mind that the values of the stochastically estimated log-likelihood are very closed from each others leading to take this result with some precaution. Finally it is clear that model M 1 does not lead to a better fit, as the log-likelihood value is only increasing minimally.

Models

Model M 2 is withheld for the current data set. The correlation structure in this model demonstrates first the substantial correlation between the udder quarter infection times within a cow. Therefore, an infected udder quarter is a higher risk factor for the three other udder quarters of the same cow than for udder quarters of a different cow. Many of the observed infections are from bacteria that are widespread in the environment so it could be rather due to the status of the cow than to the nearby presence of an infected udder quarter that makes that the other udder quarters of the same cow are more at risk.

Likelihood ratio tests

In this section we investigate through a simulation study the finite sample size property of the likelihood ratio test (LRT) for testing the hypothesis H 0 : the true model is model M 4 against the hypothesis H 1 : the true model is model M 3 . To that purpose, we repeat 200 times the following experiment: we simulate under model M 4 a data set using the general setting described in Subsection 9.4.1, we estimate from these simulated data the MLE of the parameters in model M 4 and in model M 3 and compute the corresponding log-likelihood ratio test statistic. Then we compute the empirical confidence level for different levels α ∈ {0.01, 0.05, 0.10}. Note that the model M 4 is a sub-model of model M 3 with the value of parameter ρ 1 being equal to zero which is not at the boundary of the parameter space. Thus the likelihood ratio test statistic converges in law towards a chi square random variable with one degree of freedom (see [START_REF] With Discussion | Theoretical statistics[END_REF]). The corresponding quantiles are then respectively equal to 6.63, 3.84 and 2.70. We repeat this experiment for different values of the parameter σ 2 chosen in {4, 9, 16}. The results are presented in Table 10. For example, for σ 2 = 9, we get an empirical percentage of rejection for the LRT of 0.065 which is close to the nominal level α = 0.05. We observe that the empirical percentage of rejection differs more from the nominal level α when the variance parameter equals to 16. We conduct the same experiment taking H 0 : the true model is M 3 against hypothesis H 1 : the true model is M 2 . We fix the residual variance parameter σ 2 to 9 and let the correlation parameter ρ 1 vary between 0.3, 0.5 and 0.8. The corresponding results are presented in Table 11. We observe that the empirical percentage of rejection differs more from the nominal level when the correlation parameter ρ 1 equals 0.3.

The empirical results given in Tables 10 and11 demonstrate that the LRT can be applied for a finite sample size N equal to 1196 when the parameter σ 2 is closed to 9 (and smaller) and when the parameter ρ 1 is closed to 0.5 (and larger). This heuristic observation ensures the validity of the LRT with the mastitis data in the previous section.

Conclusions and perspectives

This manuscript summarizes my contributions to the field of statistical research. They mainly deal with developments, studies and applications of new stochastic algorithmic methods derived from the Expectation Maximization algorithm for estimation in complex latent variable models. I attempted to address the theoretical aspects and the practical ones in each of my contributions. Most of my research was inspired by practical problems. Their consequences were mainly to provide more flexibility in modeling since they give solutions for parameter estimation in many complex latent variable models used in different application fields. This makes it possible, in particular, to more efficiently assess numerous complex phenomena.

All of the research topics mentioned in this document were not treated in the same way. Some were explored in depth from several different angles, leading to other developments. Other ones, generally more recent, are still being investigated and will most certainly lead to new research directions. I mentioned in the manuscript below some contributions some open questions which could lead to further developments. Besides my major actual perspectives motivated by practical questions are detailed in the next paragrafs.

Estimation in partially observed multi-type branching processes

In collaboration with Catherine Larédo (INRA, MIA), we will address the estimation in partially observed multi-type branching processes with immigration. This work is motivated by the analysis of the dynamics of rape populations. This plant has two specific properties: first, it is feral, meaning that it is able to grow by itself without particular care; second, the seeds can remain in a seed bank in the soil for several years before eventually giving rise to a new plant. A multi-type branching process is well adapted to modeling this dynamics where the different types correspond to the flowers, the seeds on the soil, the seeds in the seed bank in the soil, etc. Nevertheless, only some types were observed, namely the flowers. The model parameters such as the probability of a seed in the seek bank in the soil eventually giving rise to a new plant, characterize the dynamics and are of great interest for controlling the behavior of such a population. The estimation of these parameters is usually done in branching processes through optimization of contrast processes such as conditional least squares. However, in the presence of unobserved types in multi-type branching processes, this is no longer possible. One alternative to these approaches could be to consider the maximum likelihood approach. Note that the unobserved types are latent variable of the multi-type branching processes. However it is not possible to apply directly an EM like algorithm, the model being non exponential when considering usual reproduction and immigration distributions. Therefore we plan to consider an approximated Gaussian model for the multi-type branching process in which the estimation can be carried out through maximum likelihood by using a stochastic estimation algorithm. Our purpose is to adapt the work presented in Section 2.2 to this context to compute the MLE in the approximated Gaussian model. The proposed estimation algorithm will be apply to a rape dataset collected in Selomnes. We also plan to study the asymptotic property of the maximum likelihood estimate in the approximated model under misspecified model conditions. Our approach may be generalized to other multivariate partially observed models.

Modeling and estimation of epidemic dynamics

In collaboration with Elisabeta Vergu, (INRA, MIAJ), we are studying the epidemic dynamics modeled by compartment models and estimation in such models. The propagation of epidemic is often described by compartments, where one compartment corresponds to one possible status of an individual, e.g., sensitive, infected or recovered, in the simple case. The estimation task is difficult since observations are very often partial. Models in which the estimation can be carried out often require stringent structural assumptions of the population dynamics, e.g. Markovian. Such an assumption is of course not realistic since it implies that the recovering rate is constant, i.e., the sojourn times in the different compartments are distributed from an exponential distribution that is not realistic. We propose to consider a new non Markovian dynamics compartment model, that also takes the errors raised by the observation process into account. We plan to consider the maximum likelihood estimate and to evaluate it in practice through a stochastic estimation algorithm. We project to study the theoretical properties of the maximum likelihood estimate and of the stochastic estimation algorithm. Simulations studies and application to real epidemic data will also be carried out.

Modeling and estimation in survival data analysis

The two contributions to estimation in the frailty models presented in the third part of the manuscript have given rise to several interesting questions.

In collaboration with Charles El-Nouty (University Paris 13 Sorbonne Paris Cité, LAGA), we plan to address the case of maximum likelihood estimation in frailty models with a piecewise constant baseline function or a non parametric one. From a modeling point of view, this will make it possible to be free from the choice of a specific parametric form for the baseline, which is not obvious in practice. Another similar theoretical development will be to consider the partial likelihood, as was done for the estimation in the Cox model, rather than the marginal likelihood. This also allows us to bypass the effect of the parametric choice for the baseline. Finally, motivated by practical applications in epidemiology, we are interested in the estimation in competing risk models that integrate frailty terms. In competing risk approaches, the hazard rate is obtained as the sum of several independent hazards, each corresponding to different risks. Such modeling is, for example, well adapted to considering the risk of death, on the one hand, from a cancer, and, on the other, from cardiovascular disease.

In collaboration with Luc Duchateau and Klaartje Goethals (Faculty of Veterinary Medicine, Ghent University), we plan to assess the correlation structure in a malaria epidemic dataset in Ethiopia by introducing spatial modeling in the correlation matrix of the frailty vector. Moreover, in this dataset, the individuals are clustered through two stages, namely the villages and the houses. This will be modeled by using hierarchical models.

Another open question in the frailty model setting is the one of model validation. It would be interesting to propose goodness-of-fit tests, in particular for testing hypotheses on the covariates used in the design matrices and on the correlation structure of the frailty terms.

Analysis of plant growth

In collaboration with Hervé Monod (INRA, MIAJ), Alain Charcosset (INRA, URGV Le Moulon) and François Tardieu (INRA, LEPSE), we are interested in genotype by environment interaction. In particular, we plan to address parameter estimation in dynamic plant growth models based on several types of data including different genotypes and different environmental conditions. Indeed multiple information sources such as bibliographies, datasets obtained on platforms and datasets obtained in field conditions exist. All of this information is relevant when estimating model parameters. However, since the objective is to use the dynamic model to predict phenotypic characteristics in field conditions for existing or possibly new genotypes in different climatic scenarios. Thus parameter estimation requires close attention, particularly concerning the use of platform datasets which do not report all the phenomena existing in fields conditions. Therefore we have to integrate in the global estimation process, possibly also in the modeling, all these different types of data taking into account their specificity. Besides modeling task can also be addressed since most of the dynamic models are deterministic, although they model stochastic plants behavior. Thus it would be interesting to enrich these models by adding stochastic terms taking the existing variablity into account.

This work takes place within the project "investissement d'avenir" AMAIZING coordinated by Alain Charcosset.

In collaboration with Paul-Henry Cournède and Charlotte Baey (ECP, Digiplante), we are provided with a given deterministic dynamic model for plant growth, having parameters which may depend on the plants genotype. First, we are interested in identifying which of the model parameters could be considered fixed for several genotypes to reduce the model dimension. We plan to use non linear mixed effect models and to determine which parameters can be modeled as population parameters to analyse some arabidopsis data. In a second time, we plan to cluster the different genotypes of our dataset by using a mixture model in the spirit of [A5].

This work takes place within the project "plante entière" of the Institut de Modélisation des Sciences du Vivant (IMSV), directed by Vincent Fromion (INRA, MIG).

In collaboration with Adrienne Ressayre (INRA, URGV Le Moulon), we plan to analyze maize leave dataset. Geneticians are interested in leaf area, since it is directly link with photosynthesis and therefore with plant development, resulting in early or late flowering time. Moreover, leaf area seems to be a genotypic characteristic. To validate this hypothesis, we plan to use the BME template model presented in [A3] to estimate leaf patterns for several genotypes possibly clustered.

This work takes place within to a global project ITEMAIZE dedicated to the study of the early development of maize, in particular the floral transition, coordinated by Christine Dillmann (UPS, URGV Le Moulon) actually submitted.

Motivated by the collection processes for agronomic data in plant growth studies, I am interested in considering more precisely the particular form of partially observed data called current status data: observations are composed of the observation date and of the status of the individual at this date. This can be viewed as a form of complex censoring. Let us consider the following example: the interest is in the appearance of the third leaf of the maize plant. In practice, may be the technician crosses the field every two days only, leading to the observation, for example, that on day d, the maize plant has two leaves, and that on day d + 2, it has three leaves, generating the censored information that the third leaf appears between days d and d + 2. Such data frequently arise in plant growth studies in agronomy. Note that such data are also common in the field of epidemiology, leading to many possible applications. Modeling such partially observations leads to specific complex latent variable models. Parameter estimation is therefore particularly intricated. Existing methods are based on several types of approximations and suffer from theoretical background and numerical limitation. New subtantial developments from both the theoretical and the practical point of view are therefore required in this field.

Figure 1 .

 1 Figure 1. Ten examples of autocorrelations of the MALA (blue) and AMALA (red) samplers to target the 100 dimensional normal distribution with anisotropic covariance matrix.

Figure 2 .

 2 Figure 2. Left: twenty images per digit of the training set used for the estimation of the model parameters (inverse video). Right: same images with additive noise of variance 1.

Figure 3 .

 3 Figure 3. Estimated templates using the five algorithms on noise free and noisy data. The training set includes 20 images per digit. The dimension of the hidden variable is 72.

Figure 4 .

 4 Figure 4. Estimated templates using MALA and AMALA samplers in the stochastic EM algorithm on noisy training data. The training set includes 20 images per digit. The dimension of the hidden variable increases from 72 to 200.

Figure 5 .

 5 Figure 5. Medical image template estimation. Top rows : 10 Corpus callosum and cerebellum training images among the 47 available. Bottom row : (a) mean image. (b) FAM-EM estimated template. (c) Hybrid Gibbs -SAEM estimated template. (d) AMALA-SAEM estimated template.

Figure 6 .

 6 Figure 6. Estimated templates with varying numbers of control points (Left: 4. Middle: 9. Right: 16). Top: fixed control point model. Bottom: estimated control point model.

Figure 7 .

 7 Figure 7. Five training images from the mouse mandibles.

Figure 8 .

 8 Figure 8. Estimated templates of the mouse mandible images obtained with 260 fixed control points (left), with 117 (middle) and 70 (right) estimated control points.

Figure 9 .

 9 Figure 9. Evolution of the estimated templates and of their number of active control points with respect to the threshold parameter. From left to right: λ equals to 0.3, 0.45, 0.6, 0.75 and 0.8.

Figure 10 .

 10 Figure10. Estimated templates with their optimal numbers and positions of control points.

Figure 11 .

 11 Figure 11. Estimated prototypes of the two components model for each digit (40 images per class; 100 iterations; two components per class).

  log S i (t)) = x i β + log t 0 λ 0 (s)ds .

Figure 12 .

 12 Figure 12. Plots of t → log(-log(S(t)) in red for group A and in blue for group B.

Figure 13 .

 13 Figure 13. Plots of t → log(-log(S(t)) in red for the group without treatment and in blue for the group with treatment for a real dataset.

•

  F4 : Regularity of the baseline function. The function λ 0 belongs to the set of functions defined on R + at values in R + parameterized by the vector α taking values in an open subset A of R a , which are twice continuously differentiable on A. • F5 : Exponential family for the frailty distribution. The probability density function f η of the frailties belongs to the set of exponential probability density functions where η takes values in an open subset B of R b , which are twice continuously differentiable on B.

Theorem 8. 1 .

 1 Assume that (F1-F5), (H1-H2) and (SAEM1-SAEM3) are fulfilled. Let (θ k ) be the sequence generated by the SAEM-MCMC algorithm. We then have, with probability 1lim k→+∞ d(θ k , L) = 0,where the distance of x to the closed subset A is denoted by d(x, A) and the set of stationary points of log L obs N by L = {θ ∈ Θ, ∂ θ log L obs N (y, δ; θ) = 0}.

Figure 14 .

 14 Figure 14. Plot of t → log(-log(S(t)) in red for the group without treatment and in blue for the group with treatment.

Figure 8

 8 Figure 8.3.1 shows that the proportional hazards model assumption is not fulfilled by this dataset. We therefore consider a Gaussian frailty model with two frailty terms, referred to as Model I and defined as

Figure 15 .

 15 Figure 15. The correlation structure between the random effects of the four udder quarters for model M 1 with 'FL' front left, 'FR' front right, 'RL' rear left and 'RR' rear right udder quarter.

  Truncation step: If s is outside the current compact set K κ k-1 or too far from the previous value s k , then restart the stochastic approximation in the initial compact set, extend the truncation boundary to K κ k and start again with a bounded value of the missing variable. Otherwise, set (z k , s k ) = (z, s) and keep the truncation boundary to K κ k-1 . Update the sequence ζ k and κ k following Andrieu et al.

Table 1

 1 Comparison of EM and SAEM estimates after 100 and 1000 iterations.

	Parameters	β 1	β 2	µ	τ 2	σ 2
	θ 0	650	250	100	50	10
	θ EM 100 θ EM 1000 θ SAEM 100 θ SAEM 1000	727.89 348.06 192.05 1001.45 61.51 727.91 348.07 192.05 1001.49 61.51 725.34 346.11 191.46 1020.26 60.56 727.36 347.67 191.93 1003.76 61.52

Table 2

 2 Estimation of the standard deviation of θ EM and θ SAEM obtained after 100 and 1 000 iterations.

	Parameters	β 1	β 2	µ	τ 2	σ 2
	σ(θ EM 100 ) σ(θ EM 1000 ) σ(θ SAEM 100 σ(θ SAEM 1000	) )	13.51 13.04 14.15 633.39 14.70 13.51 13.04 14.15 633.40 14.70 12.89 12.51 13.83 604.53 13.41 13.51 13.04 14.15 633.40 14.70

Table 3

 3 Pharmacodynamic model: comparison of parameter estimates. The means and the estimated square roots of the MSEs between parentheses, based on 50 simulations.

	Parameters	Exact FOCE		LAP		EM		SAEM		σ( θ	M LE )
	µ 1	105	105.5	(1.8)	105.3	(1.6)	105.4	(1.7)	104.7	(1.5)	1.4
	µ 2	12	12.2	(1.2)	12.4	(1.3)	12.3	(1.3)	11.8	(1.3)	1.0
	µ 3	10	9.0	(2.8)	9.7	(2.4)	9.7	(1.6)	10.1	(0.9)	0.6
	τ 2 1 τ 2 2 τ 2 3	64 36 12.25	59.7 (20.9) 31.5 (11.0) 6.6 (7.6)	58.4 (20.2) 30.7 (11.9) 13.3 (6.1)	60.0 (20.1) 30.9 (10.7) 10.1 (2.9)	62.1 (16.6) 34.4 (10.8) 11.2 (3.0)	14.5 7.6 2.8

Table 4

 4 Model I. Estimation of the parameters θ = (λ 0 , β, σ 2 0 , σ 2 1 ), of the mean model-based standard error in parentheses and of the marginal log-likelihood for the bladder cancer data

	Algorithm / Estimate	λ 0 (×10 -4 )	β	σ 2 0	σ 2 1	log(L obs N )
	EM-Laplace	10.563(0.6487) -0.0167(0.0704) 0.4074(0.0923) 0.0032(0.0007) -9490.47
	MCEM	7.261(0.6441)	-0.2178(0.1224) 0.0915(0.0080) 0.1783(0.0156) -9366.59
	SAEM-MCMC	7.265(0.7186)	-0.2313(0.0880) 0.0840(0.0292) 0.1849(0.0637) -9365.41

Table 5

 5 Models I and II. Estimation of the parameters, of the mean model-based standard error in parentheses and of the marginal loglikelihood for the bladder cancer data

	Model / Estimate	λ 0 (×10 -4 )	β	σ 2 0	σ 2 1	σ 01	log(L obs N )
	I	7.265(0.7186) -0.2313(0.0880) 0.0840(0.0292) 0.1849(0.0637)	× × ××	-9365.41
	II	7.243(0.4869) -0.2540(0.0699) 0.0306(0.0002) 0.1077(0.0006) 0.0573(0.000003) -9357.51

Table 6

 6 Model I in a moderate setting with true value θ = (0.077, 0.7, 0.4, 0.8) (where λ 0 = 0.077 is known). Estimation with the three algorithms of the parameters and of standard errors in parentheses (first number: mean model based; second number: empirical)

	Algorithm / Estimate	β	σ 2 0	σ 2 1
	EM-Laplace	0.849(0.031/0.296) 0.206(0.067/0.215) 0.765(0.204/0.481)
	MCEM	0.765(0.179/0.151) 0.366(0.195/0.116) 0.778(0.289/0.195)
	SAEM-MCMC	0.729(0.072/0.122) 0.392(0.128/0.098) 0.768(0.215/0.133)

Table 7

 7 Model I in a heavy setting with true value θ = (0.077, -0.182, 0.4, 0.8) (where λ 0 = 0.077 is known). Estimation with the three algorithms of the parameters and of standard errors in parentheses (first number: mean model-based; second number: empirical)

	Algorithm / Estimate	β	σ 2 0	σ 2 1
	EM-Laplace	-0.155(0.129/0.156) 0.385(0.085/0.094) 0.766(0.165/0.200)
	MCEM	-0.167(0.189/0.174) 0.391(0.199/0.113) 0.782(0.270/0.228)
	SAEM-MCMC	-0.178(0.132/0.152) 0.392(0.120/0.109) 0.785(0.239/0.192)

Table 10

 10 Empirical percentage of rejection of the LRT of H 0 : the true model is M 4 against H 1 : the true model is M 3 at level α.Table 11Empirical percentage of rejection of the LRT of H 0 : the true model is M 3 against H 1 : the true model is M 2 at level α (σ 2 = 9).

	α/σ 2	4	9	16
	0.01	0.005 0.015 0.020
	0.05	0.045 0.065 0.080
	0.10	0.110 0.125 0.135
	α/ρ 1	0.3	0.5	0.8
	0.01	0.030 0.015 0.010
	0.05	0.085 0.060 0.055
	0.10	0.135 0.115 0.095

Scientific path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Scientific context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

follows:

λ ij (t|b i ) = λ 0 (t) exp b 0i + x t ij (β + b 1i ) ,

for 1 ≤ i ≤ N, 1 ≤ j ≤ n i and t ≥ 0, where b i = (b 0i , b 1i ) t , the random variables (b 0i ) 1≤i≤N (respectively, (b 1i ) 1≤i≤N ) are independent and identically distributed Gaussian N (0, σ 2 0 ) (respectively, N (0, σ 2 1 )), and the sequences (b 0i ) 1≤i≤N and (b 1i ) 1≤i≤N are independent. Note that b 0i and b 1i do not have a symmetric role. Model I, already considered by [START_REF] Abrahantes | Comparison of different estimation procedures for proportional hazards model with random effects[END_REF] and by [START_REF] Legrand | A Bayesian approach to jointly estimate centre and treatment by centre heterogeneity in a proportional hazards model[END_REF], is used to analyze the dataset of bladder cancer and also to carry out our simulation studies in the next section. The event times T ij correspond to times from randomization to the date of the first bladder recurrence, censoring the patients without recurrence at the date of the last available follow-up cystoscopy. The covariate x ij is equal to zero if the j-th patient of the i-th group receives no further intra-vesical treatment and equal to one otherwise. More generally, this covariate indicates whether the patient is classified in the good or poor prognosis group based on the particular prognostic index considered (see [START_REF] Legrand | Validation of prognostic indices using the frailty model[END_REF]). The variable b 0i can be understood as a random center effect, whereas the variable b 1i is viewed as a random treatment per center interaction. Finally, the parameter β is the fixed treatment effect. We refer to [START_REF] Sylvester | Predicting recurrence and progression in individual patients with stage ta t1 bladder cancer using eortc risk tables: A combined analysis of 2596 patients from seven eortc trials[END_REF] for detailed and complete information on the dataset.

Before starting the numerical studies, we have to choose the parametric baseline function λ 0 . The SAEM-MCMC algorithm can handle a large class of functions. For example we may choose the Gompertz function λ 0 (t) = λ exp(γt), λ > 0, γ ∈ R, or the Weibull function λ 0 (t) = λ ρ t ρ-1 , λ > 0, ρ > 0 for t ≥ 0 (see Duchateau and Janssen (2008, p 29)). Nevertheless, as in [START_REF] Abrahantes | Comparison of different estimation procedures for proportional hazards model with random effects[END_REF], we assume in this section that the baseline function is constant to compare our results with others found in the literature. Hence, the vector of unknown parameters becomes: θ = (λ 0 , β, σ 2 0 , σ 2 1 ).

We can easily verify that the regularity assumptions (F1-F5), (H1-H2) and (SAEM2) are fulfilled by Model I. To satisfy the assumption (SAEM1), we choose the sequence (γ k ) k so that if 1 ≤ k ≤ 50, then γ k = 1; otherwise: γ k = 1 (k-50) 2/3 . Finally, we choose the transition probability Π θ as a hybrid Gibbs sampler, also known as the Metropolis Hastings within Gibbs algorithm, in order to verify assumption (SAEM3). The proposal distribution is chosen to be equal to the distribution of the frailty terms.

We apply three algorithms to estimate the parameters θ = (λ 0 , β, σ 2 0 , σ 2 1 ). We emphasize that we implement the MCEM algorithm (respectively, the EM-Laplace algorithm) according to the method given in [START_REF] Ripatti | Maximum likelihood inference for multivariate frailty models using an automated Monte Carlo EM algorithm[END_REF] (respectively, in Cortiñas [START_REF] Cortiñas Abrahantes | A version of the EM algorithm for proportional hazard model with random effects[END_REF]).

Different random initializations were tested. To investigate the nature of the limit points, we compute an estimate of the logarithm of the marginal likelihood for each limit point using a Monte Carlo sum. We present the relevant numerical results in Table 4 corresponding to the trajectory that gives the estimate of the logarithm of the marginal likelihood with the biggest value.

Let us first comment on the behavior of each algorithm, on the one hand, with respect to the convergence of the trajectories and, on the other, with respect to the computing time.

Consider the convergence problem first. When we apply the SAEM-MCMC algorithm or the MCEM algorithm, all the trajectories converge. The situation is totally different with the EM-Laplace algorithm.

The matrix Σ being positive definite, it follows:

We define the sub-model M 4 of model M 1 by assuming that there is no correlation between any pair of random effects. Note that the model M 4 is a sub-model of model M 3 obtained by setting the value of ρ 1 equal to zero.

The covariance matrix Σ of model M 4 is given by:

where σ 2 > 0 and I 4 denotes the identity matrix of size 4. Model M 4 thus assumes independence between the event times in a cluster, or stated differently: the fact that an udder quarter has an infection does not have an influence on the hazard of infection for the three other udder quarters. Model M 4 is termed as univariate frailty model, because each subject has its own frailty term.

We will now focus on parameter estimation in these four frailty models. We insist again on the fact that the four models M 1 , M 2 , M 3 and M 4 are nested. This allows us to apply the likelihood ratio test to compare them in the forthcoming sections.

Statistical analysis of the mastitis dataset

We choose a frequentist approach and consider the Maximum Likelihood Estimate (MLE) for the parameter θ, namely the value θN of θ which maximizes the marginal likelihood L obs N (Y, ∆; θ) (see [START_REF] Duchateau | The Frailty Model[END_REF]). We apply the estimation algorithm for frailty models presented in [A6].

We focus now on the theoretical properties of the algorithm. Some particular assumptions on the frailty model have been proposed and stated in Theorem 1 in [A6] to ensure the almost sure convergence of this stochastic estimation algorithm toward the set of stationary points of the observed likelihood. They are fulfilled by the Weibull baseline function and by the multivariate Gaussian distribution for the frailty. Moreover the model M 1 , as well as the three sub-models M 2 , M 3 and M 4 , satisfy the other additional regularity assumptions (denoted by (H1-H2) and (SAEM2) in Theorem 1) required to apply the convergence theorem. Furthermore we choose a step size sequence ν = (ν k ) k and a transition probability kernel Π θ for the stochastic approximation procedure which satisfy assumptions (SAEM1) and (SAEM3) of Theorem 1. Thus applying the result of Theorem 1 in [A6] proves that the parameter estimated sequence (θ k ) k generated by our algorithm converges almost surely towards the set of stationary points of the observed likelihood Parameters and mean model based standard errors estimates are obtained applying the methods described in Section 8, for each of the four models (cf Table 8). The mean model based standard errors are obtained as the square root of the diagonal components of the inverse of the Fisher Information Matrix.

Table 9

The mean of the parameter estimates of the 500 generated datasets for the different parameters with the empirical and mean model based standard errors in parenthesis in models M 2 and M 3 . On the other hand, model M 2 also reveals that the risk factor for an udder quarter is increased more if the udder quarter is in the same region, i.e., either front or rear, as the infected udder quarter. This finding means that the nearby presence of an infected udder quarter has the effect of increasing the risk of infection.

Simulation studies

The aim of the simulation studies is the validation of the results obtained for the mastitis data set in Section 9.1. In the first subsection 9.4.1 the bias of the parameter estimates for the two best models M 2 and M 3 is investigated through simulations in similar settings as the real dataset. In the second subsection 9.4.2, the finite sample size property of the likelihood ratio test for models M 2 , M 3 and M 4 is considered to ensure that the level of the test is close to the nominal one. 9.4.1. Parameter Estimation Evaluation for Models M 2 and M 3 Data are simulated from models M 2 and M 3 using a similar setting as in the mastitis dataset. In total N = 1196 clusters are chosen, each cluster consisting in 4 infection times. Also the effect of parity and location is included.

The following parameter values are used to generate the data: λ 0 = 0.01, γ = 2.75, β 1 = 0.90, β 2 = -0.32, σ 2 = 9, ρ 1 = 0.8 and ρ 2 = 0.7 (only for model M 2 ).

For each of the two models M 2 and M 3 , 500 datasets were generated and the methods of Section 8 are applied to obtain the parameter and mean model based standard errors estimates for each of the datasets. The simulation results are presented in Table 9 and summarized by the mean of the parameter estimates of the 500 datasets, the empirical standar errors and the mean of the mean model based standard errors of the 500 datasets.

The mean of the estimates is closed to the true value for all parameters in each of the two models M 2 and M 3 . In most cases, the values obtained for the estimated variances are also coherent, the model based one being lower than the empirical one.

Bibliography

List of works

List of publications