# Analyse longitudinale multivariée par modèles mixtes et application à l'épidémie de la malaria

Eric Houngla Adjakossa

► To cite this version:

Eric Houngla Adjakossa. Analyse longitudinale multivariée par modèles mixtes et application à l'épidémie de la malaria. Mathématiques [math]. Université Pierre et Marie Curie, 2017. Français.
NNT: 2017PA066014 . tel-04021346

HAL Id: tel-04021346

https://theses.hal.science/tel-04021346

Submitted on 9 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THÈSE

présentée pour obtenir

## LE GRADE DE DOCTEUR EN SCIENCES DES UNIVERSITÉS PIERRE ET MARIE CURIE ET D'ABOMEY-CALAVI

Spécialité : Mathématiques Appliquées

par

Eric Houngla ADJAKOSSA

## Analyse longitudinale multivariée par modèles mixtes et application à l'épidémie de la malaria

dirigée par MM. Grégory NUEL et M. Norbert HOUNKONNOU

Rapporteurs : M. Daniel **COMMENGES** Université de Bordeaux
M. Romain L. **GLELE KAKAÏ** Université d'Abomey-Calavi

Soutenue le 03 avril 2017 à l'Université Pierre et Marie Curie devant le jury composé de

<table><tr><td>M. Daniel</td><td><strong>COMMENGES</strong></td><td>Université de Bordeaux</td><td>Président</td></tr><tr><td>M. Amaury</td><td><strong>LAMBERT</strong></td><td>Université Pierre et Marie Curie</td><td>Examineur</td></tr><tr><td>M. Grégory</td><td><strong>NUEL</strong></td><td>Université Pierre et Marie Curie</td><td>Directeur</td></tr><tr><td>M. M. Norbert</td><td><strong>HOUNKONNOU</strong></td><td>Université d'Abomey-Calavi</td><td>Directeur</td></tr><tr><td>M. Romain L.</td><td><strong>GLELE KAKAÏ</strong></td><td>Université d'Abomey-Calavi</td><td>Rapporteur</td></tr><tr><td>M. David</td><td><strong>COURTIN</strong></td><td>Université Paris Descartes</td><td>Examineur</td></tr></table>
"The only thing we have to fear is fear itself."

Franklin D. Roosevelt

Je dédie cette thèse à tous ceux qui croient en l'AMOUR
Remerciements

Je tiens, en premier lieu, à exprimer toute ma gratitude à mes deux directeurs de thèse Grégory NUEL et Mahouton Norbert HOUNKONNOU pour avoir accepté de m’accompagner et de me soutenir pendant les trois années qu’ont duré mes travaux de thèse. Lors de nos multiples rendez-vous, j’ai pu profiter de vos nombreuses capacités scientifiques. Un merci tout particulier à Grégory avec qui j’ai passé plus de temps pendant toutes ces années, et qui a su me rendre progressivement indépendant, en tout cas beaucoup plus que je ne l’étais en tout début de thèse. Très sincèrement merci.

Merci à tous les membres du jury pour l’honneur qu’ils me font d’être présents aujourd’hui. Je suis très reconnaissant envers M. Daniel COMMENGES et M. Romain Lucas GLÈLÈ KAKAÏ AGBIDINOUKOUN d’avoir pris le temps de lire de façon critique mon manuscrit de thèse en tant que rapporteurs, et ce malgré leur emploi du temps déjà bien chargé et leurs obligations personnelles.

Mes remerciements vont également à toute l’équipe de l’UMR 216 de l’IRD qui a généré et mis à notre disposition les données d’anticorps anti-palustres sur lesquelles j’ai travaillé. Remerciements particuliers à André GARCIA et à Gilles COTTRELL pour leur soutien depuis mon Master 2. Je vous en suis profondément reconnaissant. Un merci particulier à David COURTIN qui, depuis quelques années, dès que j’en exprime le besoin, trouve toujours le temps de m’expliquer très clairement différentes notions biologiques. Merci aussi d’avoir accepté d’être présent dans le jury, toujours pour m’accompagner jusqu’au bout. J’en suis sincèrement touché.

Toute ma gratitude va également à Amaury LAMBERT pour l’honneur qu’il me fait d’être membre de mon jury de thèse malgré son emploi du temps très chargé. Sincèrement merci.

Je tiens à exprimer ma gratitude aux Professeurs intervenus dans le Master STAFAV (Statistique pour l'Afrique Francophone et Application au Vivant) de Cotonou dont je suis un pur produit. Gratitude particulière aux Professeurs Jean-Marc BARDET, Olivier WINTENBERGER et Simplec DOSSOU-GBÉTÉ. Mes remerciements vont également aux autorités ainsi qu’à l’administration de la CIPMA chaire-UNESCO à Cotonou pour m’avoir beaucoup aidé dans les diverses démarches administratives. Je pense au Président Norbert HOUNKONNOU, au Dr BALOÏTCHA ainsi qu’aux secrétaires.

Je remercie le LPMA (Laboratoire de Probabilités et Modèles Aléatoires), mon laboratoire d’accueil à Paris, où j’ai pu bénéficier de bonnes conditions de travail. Je pense à Florence Deschamps, Josette SAMAN, Serena BENASSU, Altaïr PELISSIER et Khashayar DADRAS. Je remercie également tout le personnel de l’Ecole Doctorale de Sciences Mathématiques de Paris Centre 386, et plus particulièrement Jean-François VENUTI et Corentin LACOMBE.
Pendant ces années de thèse, les doctorants du LPMA ont joué un grand rôle dans ma vie et ont été pour moi une seconde famille. Un grand merci à mon "club" de footing composé de: Pierre Antoine, Alice, Sarah, Paul, Carlo et Nicolas le chef d'équipe. Je remercie également mes collègues de bureau avec qui j'ai passé d'agréables moments. Un merci tout particulier à Olga et Omar dont la silencieuse sympathie continue de m'épater. Un grand merci à Nelo qui n'a cessé de m'écouter patiemment, de m'encourager et de m'aider au besoin. Des moments comme le ciné club des thésards sont à jamais gravés en moi. Du fond du cœur, je vous exprime, à chacun et à tous, toute ma gratitude.

Je tiens à remercier tous les enseignants que j'ai connus depuis que j'ai mis pieds à l'école. Quelques uns parmi eux sont définitivement restés gravés dans ma mémoire. Merci de tout cœur à Barnabé HOUNGBEME qui a été mon maître en classe de CM2, et dont la rigueur et la pédagogie m'ont beaucoup apporté et ont semé en moi l'envie de réussir. Profonde gratitude aussi à Gustave GODONOU, mon professeur de français en classe de 3e, dont l'acharnement à corriger la qualité de notre rédaction est aussi resté gravé en moi. Je me rappelle encore toutes ses stratégies pour nous aider; tellement cela m'a marqué ! Mes sincères remerciements à Cyrille HOUNKPATIN, mon professeur de mathématiques en classe de terminale, dont l'humilité et la pédagogie m'ont définitivement rendu confiant.

Mes remerciements vont aussi à toute la communauté des ECKistes, béninois comme français. Leur amour m'a constamment soutenu et nourri. Merci à Thérèse, Achille, Clovis, Estelle, Chimène, Bernadette, Dovis, Sonia, Marie-Agnès et Doune. Merci aussi aux familles AGLI, HOUNWANOU, KOUTON, HOUNGBO, NOUHOUN, de SOUZA (de Hugor) et ADEOCHUN (de Djahafar). Un merci tout particulier à Sabira EZIN que je ne finirai jamais de remercier pour son AMOUR. Profonde gratitude aussi à tous les Maîtres ECK qui m'ont comblé de leur soutien et de leur AMOUR. Merci de tout cœur.

Mon envie de faire de la recherche scientifique aurait pu être étouffée sans le secours de personnes au grand cœur comme Jacques SEGLA et Luc HESSOU. Je vous suis infiniment reconnaissant pour votre soutien déterminant sans lequel je ne serais peut-être pas en train de finir cette thèse en ce moment. Pour moi, vous êtes des exemples à suivre pour bâtir un Bénin plus grand, voire une Afrique plus grande.

Cette thèse, évidemment, n'aurait été possible sans l'assistance de mes parents. Merci à toi papa pour ta rigueur et ta discipline absolues sans lesquelles je ne serais pas là. Tu savais bien de quoi tu parlais quand tu me disais il y a quelques années: "Eric, tu ne peux plus jamais te laisser aller à la paresse". Merci à toi maman pour ta tendresse et tes nombreux sacrifices pour tes enfants. Puisse le temps nous donner la force de te témoigner concrètement et suffisamment notre gratitude. Merci à tous mes frères et sœurs. Merci à tous mes amis que je ne peux distinctement citer ici.

Enfin, il n'aurait pas été pensable que je vous oublie. Merci à Sandra, Shariy et Huson pour votre amour, votre soutien, votre patience et cette amitié totale. Merci pour ce qu'il y a de meilleur: l'AMOUR ici et maintenant !
# Préambule

## Contexte et objectifs

La première trace de paludisme est la présence d'ADN de *Plasmodium falciparum*⁽¹⁾ chez des momies datées de 3200 avant Jésus Christ (JC) [118]. Historiquement, l'origine du paludisme chez l'Homme serait due à une transmission venant du gorille, excluant ainsi une transmission due aux chimpanzés ou aux ancêtres de l'Homme lui-même [127]. De nombreuses références au paludisme existent dans les écrits védiques⁽²⁾ (1600 avant JC) et sur les tablettes d'argile d'Ashurbanipal⁽³⁾ (669 avant JC) en Mésopotamie [43]. Le terme "malaria" venu de l'italien (ancien) est une combinaison des mots *mala* signifiant "mauvais, insalubre" et "aria" signifiant "air", et est utilisé comme nom de maladie depuis le 17ᵉ siècle. D'ailleurs, Hippocrate, vraisemblablement parmi les premiers à décrire certains aspects et symptômes de la maladie (fièvre, aspect saisonnier), a établi une relation avec les eaux stagnantes des marais [130]. Au 19ᵉ siècle, le terme "paludisme" apparait venant du mot *palud*, lui-même dérivé du latin *palus* signifiant "marais".

Des travaux du français Charles Laveran en 1880 à ceux du britanique Ronald Ross en 1897, il a été prouvé que le paludisme est une maladie infectieuse à transmission vectorielle faisant intervenir trois acteurs. L'homme, jouant le rôle d'hôte, est infecté par un protozoaire parasite du genre *Plasmodium* qui lui a été transmis par la piqûre d'un moustique (femelle) du genre *Anopheles* [29]. La maladie se caractérise par des épisodes fébriles aigus et peut être mortelle. Des symptômes tels que maux de tête, fièvre, vomissements et frissons apparaissent au bout de sept jours ou plus après la piqûre infectante du moustique. L'intensité de la transmission du paludisme dépend de facteurs liés au parasite, au vecteur (moustique), à l'hôte humain et à l'environnement [36, 136, 176].

Les efforts mis en œuvre pour la maîtrise du développement du paludisme ont réduit sa distribution géographique mondiale de 53% en 1900 à 27% en 2002 [89]. En effet, le paludisme a commencé à régresser en Europe dès le 18ᵉ siècle, puis en Amérique du Nord au milieu du 19ᵉ siècle. Il a disparu en Angleterre, en France ainsi que dans de nombreux pays européens vers la fin du 19ᵉ siècle. Il se maintenait par contre en Italie, en Grèce, en Allemagne et sévissait surtout dans les régions tropicales et subtropicales [130].

Malgré les efforts acharnés fournis depuis des décennies, le paludisme est resté un problème d'ampleur mondiale qui se concentre sur la santé et le potentiel économique des communautés les

<footnote>⁽¹⁾ *Plasmodium falciparum* est une des espèces de *Plasmodium*, parasites qui causent le paludisme chez l'être humain. Il est le plus dangereux de ces parasites causant le paludisme car il entraîne le taux de mortalité le plus élevé.</footnote>

<footnote>⁽²⁾ Le Veda, du sanskrit "vision ou connaissance", est un ensemble de textes qui, selon la tradition, ont été révélés par l'audition aux sages indiens nommés *Richi* [134].</footnote>

<footnote>⁽³⁾ Tablettes réunies en faisceau par des bandelettes et portant des inscriptions à caractère prophétique [128].</footnote>
plus pauvres de la planète. C'est pour cela qu'en 2000, la cible 6C des Objectifs du Millénaire pour le Développement (OMD) liée au paludisme appelait à avoir maîtrisé cette maladie en 2015 et commencé à inverser la tendance (de 2000). Malheureusement, en 2012, selon les estimations de l'Organisation Mondiale de la Santé (OMS), 207 millions de cas sont survenus et la maladie a tué quelques 627.000 personnes - pour la plupart, des enfants de moins de cinq ans vivant en Afrique. En moyenne, le paludisme tue un enfant chaque minute [138]. Le Rapport sur le paludisme dans le monde 2015 résume les progrès accomplis en matière de contrôle et d'élimination de la maladie à la date-butoir de 2015 et mentionne que "malgré des progrès remarquables, il reste beaucoup à faire". La Stratégie technique mondiale de lutte contre le paludisme 2016-2030, approuvée par l'Assemblée mondiale de la Santé en mai 2015, définit des objectifs ambitieux pour 2030, notamment réduire d'au moins 90% l'incidence du paludisme et la mortalité associée. La version complète du Rapport sur le paludisme dans le monde 2015 est disponible en anglais à l'adresse http://www.who.int/malaria/publications/world-malaria-report-2015/report/en/.

En Afrique où surviennent le plus de cas, l'objectif est de contrôler le paludisme par des approches intégrées combinant le diagnostic, un traitement efficace et précoce, la lutte anti-vectorielle et la chimioprophylaxie, notamment par le Traitement Préventif Intermittent (TPI) des enfants et des femmes enceintes [43]. Un vaccin antipaludique serait bien entendu un apport essentiel dans cette lutte [157].

Les deux populations majoritairement à risque restent les enfants de moins de cinq ans et les femmes enceintes. C'est pour cette raison que l'Institut de Recherche pour le Développement (IRD) et en particulier l'unité UMR 216 "Mère et enfant face aux infections tropicales" s'intéressent à ces populations à travers plusieurs suivis de cohortes en zone subsaharienne (Bénin, Sénégal) dont l'objectif est de comprendre les déterminants immunologiques, parasitologiques, génétiques et environnementaux de la survenue des infections palustres dans le but d'élaborer à terme un vaccin contre le paludisme. Les enquêtes épidémiologiques menées sur ces cohortes génèrent d'importants jeux de données comportant plusieurs dizaines de variables épidémiologiques, biologiques, entomologiques, environnementaux, etc. Ces données sont le plus souvent répétées dans le temps (longitudinales) et présentent des structures de corrélation particulières dans le temps et/ou dans l'espace qu'il faut prendre en compte dans les analyses au travers d'ajustement de modèles statistiques adaptés.

Spécifiquement, les données dont l'analyse est à l'origine de cette thèse proviennent des travaux d'un projet de l'UMR 216 dont l'objectif principal était d'étudier la construction de la réponse immunitaire humorale dirigée contre *P. falciparum* durant les dix huit premiers mois de vie chez 600 enfants vivant dans une zone rurale d'endémie palustre à Tori-Bossito dans le Sud-Ouest du Bénin. Ces données incluent sept antigènes des stades sanguins asexués de *P. falciparum*, retenus en leur qualité de candidats vaccins: MSP1, MSP2 (3D7 et FC27), MSP3, AMA1, GLURP (R0 et R2). Ces antigènes interviennent de façon déterminante dans le développement d'une immunité protectrice dans la mesure où ils induisent la production d'anticorps associés à la protection clinique, notamment les immunoglobulines (Ig)G de type cytophile IgG1 et IgG3 [43].

Dans les modèles statistiques utilisés classiquement pour analyser ces données en prenant en compte leur structure, chaque profile d'anticorps était en général expliqué par quelques uns ou tous les autres profiles d'anticorps et/ou en présence d'autres facteurs (environnement, anémie, etc.).
Plusieurs modèles indépendants étaient alors ajustés à ces données en vue d'apporter des réponses adéquates aux différentes questions. Or tous ces antigènes inter-agissent au même moment dans l'organisme de l'enfant et il semble plus approprié de procéder à l'analyse de la distribution jointe de ces profiles d'anticorps, malgré l'utilité incontestable des analyses séparées (facilité d'usage, résultats aisément interprétables).

Les objectifs de cette thèse étaient donc d'apporter un support en statistique ainsi qu'une profonde connaissance de la littérature dans le but d'améliorer la compréhension des méthodes d'analyses de données longitudinales multidimensionnelles. Ils étaient également de contribuer au passage aux développements méthodologiques en vue d'une analyse plus efficace de données multidimensionnelles multi-niveaux comme celles issues des suivis épidémiologiques de cohortes relatives à la lutte contre le paludisme.

Cette thèse a été effectuée dans un contexte de co-tutelle internationale entre le Bénin et la France, précisément les Universités Pierre et Marie Curie (UPMC) et d'Abomey-Calavi (UAC). Les séjours en France aussi bien qu'au Bénin ont été conjointement financés par le Service de Coopération et d'Action Culturelle (SCAC) de l'ambassade de France au Bénin et l'Institut de Recherche pour le Développement (IRD). L'ensemble de ce travail a été réalisé sous l'encadrement scientifique des Professeurs Grégory NUEL (LPMA/UPMC) et M. Norbert HOUNKONNOU (CIPMA/UAC).

## Plan de la thèse

Cette thèse porte sur l'estimation et la sélection dans le modèle linéaire à effets mixtes et s'organise autour de cinq chapitres indépendants. Chaque chapitre s'ouvre sur un énoncé introductif des différents points qui y sont abordés et se conclut sur une synthèse.

Le chapitre 1 est une introduction générale qui présente le développement historique des méthodes d'analyse de données longitudinales et définit le modèle statistique linéaire à effets mixtes à partir d'un exemple simple. Il présente ensuite les principales méthodes d'estimation des paramètres du modèle, précise les différentes contributions scientifiques faites dans le cadre de cette thèse suivies des perspectives envisagées.

Le chapitre 2 présente un article de recherche publié dans PLOS ONE portant sur l'estimation des paramètres du modèle linéaire multidimensionnel à effets mixtes, utilisant l'algorithme EM. Les estimateurs EM présentés dans ce chapitre ont des expressions plus générales que celles rencontrées dans la littérature traitant des données longitudinales multi-variées. Ces estimateurs tels que présentés intègrent le cadre plus général d'analyse de données multidimensionnelles multi-niveaux. Dans ce chapitre, avons-nous également introduit un test de corrélation bi-varié permettant de tester la significativité de l'ensemble des corrélations entre les effets aléatoires de deux dimensions du modèle. Ce qui permettra de tester si deux composantes du vecteur de réponses (les variables dépendantes) sont à analyser séparément ou non. Ce test de corrélation bi-varié permettra également de construire un modèle multidimensionnel plus parcimonieux en terme de nombre de composantes de variance des effets aléatoires, à travers une procédure de sélection pas-à-pas ascendante.
Le chapitre 3 présente un autre article de recherche que nous avons soumis à une revue sci-
entifique. Le travail exposé ici est une généralisation de la méthode d'estimation utilisée dans le
package lme4 du logiciel R, proposé par Douglas Bates et ses collaborateurs. Au travers d'une
'étude de simulation, nous avons montré que la procédure proposée est beaucoup moins sensible
au point de départ et plus rapide, comparée à l'estimation par l'algorithme EM. Aussi fournit-elle
des estimations consistantes sur nos données simulées.

Le chapitre 4 présente une procédure de sélection d'effets fixes dans le modèle linéaire à effets mixtes. Cette procédure de sélection que nous avons appelée l'*adaptive ridge* itérative utilise une pénalité de type *L*₂ sur le vecteur des effets fixes en présence d'une matrice de poids itérativement mise à jour, et où nous utilisons la log-vraisemblance profilée des composantes de variance des effets aléatoires et des effets fixes, conditionnellement aux observations. Cette approche permet d'approcher les performances en matière de sélection d'une pénalité de type *L*₀.

Dans le chapitre 5, nous présentons davantage le cadre d'application qu'est le notre dans cette
thèse avec les données d'anticorps anti-paludisme. Ici, les différentes illustrations de nos méthodes
sur ces données sont également exposées.
# Table des matières

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1</td><td>Développement historique des modèles linéaires pour l'analyse de données longitudinales</td><td>2</td></tr><tr><td>1.2</td><td>Modèle linéaire à effets mixtes : différentes méthodes d'estimation des paramètres</td><td>5</td></tr><tr><td>1.2.1</td><td>Définition du modèle</td><td>5</td></tr><tr><td>1.2.2</td><td>Estimation des paramètres du modèle linéaire à effets mixtes</td><td>8</td></tr><tr><td>1.3</td><td>Contributions méthodologiques de la thèse</td><td>17</td></tr><tr><td>2</td><td>Multivariate Longitudinal Analysis with Bivariate Correlation Test</td><td>27</td></tr><tr><td>2.1</td><td>Introduction</td><td>28</td></tr><tr><td>2.2</td><td>Materials and Methods</td><td>30</td></tr><tr><td>2.2.1</td><td>Previous works</td><td>30</td></tr><tr><td>2.2.2</td><td>Model and notations</td><td>34</td></tr><tr><td>2.2.3</td><td>EM estimation</td><td>38</td></tr><tr><td>2.2.4</td><td>Test of the significance of Cor̂(γ₁, γ₂)</td><td>40</td></tr><tr><td>2.3</td><td>Results and Discussion</td><td>41</td></tr><tr><td>2.3.1</td><td>Simulation studies</td><td>41</td></tr><tr><td>2.3.2</td><td>Applications on real data sets</td><td>48</td></tr><tr><td>2.4</td><td>Conclusion</td><td>54</td></tr><tr><td>3</td><td>Profiled deviance for the multivariate linear mixed-effects model fitting</td><td>59</td></tr><tr><td>3.1</td><td>Introduction</td><td>60</td></tr><tr><td>3.2</td><td>Multivariate linear mixed-effects model</td><td>61</td></tr><tr><td>3.3</td><td>Parameters' estimates</td><td>63</td></tr><tr><td>3.3.1</td><td>ML criterion</td><td>63</td></tr><tr><td>3.3.2</td><td>REML criterion</td><td>67</td></tr><tr><td>3.4</td><td>Simulation studies</td><td>68</td></tr><tr><td>3.4.1</td><td>Estimates' performances</td><td>69</td></tr><tr><td>3.4.2</td><td>Comparison with EM-based estimates</td><td>71</td></tr><tr><td>3.5</td><td>Application on malaria dataset</td><td>74</td></tr><tr><td>3.5.1</td><td>Data description</td><td>74</td></tr><tr><td>3.5.2</td><td>Data analysis</td><td>75</td></tr><tr><td>3.6</td><td>Conclusion</td><td>76</td></tr></tbody></table>
INTRODUCTION

Il est courant que les observations ou les mesures faites dans plusieurs domaines d'application
des mathématiques (biologie, sociologie, statistique, économie, etc.) fassent émerger une certaine
structure de corrélation dont la compréhension et la prise en compte sont nécessaires pour l'ana-
lyse de l'ensemble des données engendrées. Ainsi parle-t-on de données corrélées. Les termes
génériques données corrélées embracent une multitude de structures de données telles que : les
observations multivariées, les données groupées, les mesures répétées, les données de type longi-
tudinal et les données spatialement corrélées [193]. Dans cette thèse, nous nous focalisons sur
l'analyse statistique de données de type longitudinal multidimensionnel qu'on appelle parfois
données de type longitudinal multivarié. Ces données proviennent des études longitudinales dont
la caractéristique essentielle est qu'un ou plusieurs attributs descripteurs des sujets (ou unités sta-
tistiques) à l'étude sont mesurés de façon répétée au cours du temps. Par exemple, on pourrait
s'intéresser à l'évolution du poids et de la taille d'enfants nés dans un hôpital un jour donné, en les
mesurant tous les trois mois pendant cinq ans. Ceci constitue un exemple de données longitudi-
nales, et si les mesures (poids ou tailles) provenant d'un même enfant sont généralement corrélées,
celles provenant de deux différents enfants peuvent ne pas l'être. Les données de type longitudinal
peuvent être recueillies de façon prospective, les unités statistiques étant suivies progressivement
dans le temps, ou de façon rétrospective, en extrayant plusieurs informations des archives rela-
tives à chaque unité statistique. Les études longitudinales, contrairement aux études transversales
où chaque unité statistique est observée (ou mesurée) une seule fois, ont l'avantage d'aider à dis-
tinguer les changements parmi les unités statistiques au fil du temps (effets du vieillissement par
exemple) [45]. Dans les études longitudinales, bien que certaines questions puissent être résolues
en analysant séparément différentes variables d'intérêt, d'autres quant à elles ne peuvent l'être
qu'en analysant conjointement ces variables [194]. Ce qui revient à considérer un vecteur dont
les composantes sont les différentes variables d'intérêt de type longitudinal : d'où la terminolo-
gie données longitudinales multidimensionnelles ou multivariées. Nous nous intéressons ici au cas
où ces variables sont des nombres réels et aux méthodes statistiques permettant de les analyser.
Parmi les modèles statistiques permettant d'analyser efficacement les données longitudinales mul-
tivariées, nous nous intéressons plus précisément au modèle linéaire à effets mixtes. La raison de
cet intérêt est expliquée à travers le développement historique des modèles linéaires pour l'analyse
de données longitudinales.
## 1.1 Développement historique des modèles linéaires pour l’analyse de données longitudinales

Un des plus anciens exemples d’analyse de données longitudinales fourni par la littérature nous vient des travaux de l’astrologue George Biddel Airy [2, Partie IV] qui portaient sur la modélisation de plusieurs niveaux d’erreurs contenues dans les mesures successives d’une même quantité supposée invariante, comme par exemple le diamètre angulaire apparent d’une planète ou encore la distance angulaire entre deux étoiles. Plusieurs mesures étaient faites quotidiennement et ce pendant plusieurs jours. Airy constata que même après avoir appliqué tous les types de correction d’erreur de mesure connus, il apparaissait parfois que les résultats obtenus différaient grandement d’un jour à l’autre. C’est pourquoi il émit l’hypothèse de l’existence de deux types d’erreur suivant le modèle

$$y_{ij} = \mu + \gamma_i + \varepsilon_{ij}, \qquad (1.1)$$

où $y_{ij}$ est la $j^e$ mesure du $i^e$ jour et $\mu$ la moyenne générale ou la “vraie” valeur de l’objet mesuré. La quantité $\gamma_i$ représente, quant à elle, l’erreur (constante) commise sur toutes les mesures effectuées le $i^e$ jour, et $\varepsilon_{ij}$ l’erreur due aux conditions particulières d’observation (atmosphère, instruments de mesure, etc.) dans lesquelles la $j^e$ mesure du $i^e$ jour a été effectuée. Plus tard, Henry Scheffé [162], proposant une alternative pour le modèle d’analyse de la variance tel qu’introduit par R. A. Fisher [60], dans le but d’analyser des données issues de dispositifs en blocs randomisés, introduisit les termes d’effet fixe et d’effet aléatoire. De façon plus explicite, selon Scheffé, le modèle d’erreurs d’Airy (équation (1.1)) s’écrit :

$$y_{ij} = \mu + \gamma_i + \varepsilon_{ij}, \text{ avec } \gamma_i \sim \mathcal{N}(0, \sigma_{\gamma}^2), \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma_{\varepsilon}^2) \text{ et } \gamma_i \perp \varepsilon_{ij}, \qquad (1.2)$$

où $\gamma_i$ s’interprète désormais comme étant l’effet aléatoire du $i^e$ jour sur la valeur mesurée $y_{ij}$, et $\varepsilon_{ij}$ le terme d’erreur (ou résiduel) du modèle. Le terme $\gamma_i$ est appelé effet parce qu’il représente une déviation par rapport à la moyenne générale $\mu$, la moyenne des mesures effectuées le $i^e$ jour étant $\mu + \gamma_i$. Les jours pendant lesquels les mesures ont été effectuées peuvent être considérés comme un échantillon aléatoire de l’ensemble $J$ des jours pour lesquels on voudrait inférer des informations. Par exemple, $J$ pourrait être l’ensemble des jours depuis la naissance du système solaire jusqu’à son extinction. En supposant que $J$ est infini et que l’estimation et l’inférence sont faites sur $J$ et non uniquement sur les jours où des mesures ont été effectuées, l’effet du jour $i$ est aléatoire (par exemple, $\gamma_i \sim \mathcal{N}(0, \sigma_{\gamma}^2)$) et non fixe. Henry Scheffé introduisit ainsi les concepts de modèle à effets fixes, de modèle à effets aléatoires et de modèle à effets mixtes [162] pour l’analyse de données longitudinales et plus généralement, l’analyse de données à structure hiérarchique, induisant une structure de correlation plus ou moins complexe à prendre en compte. Le modèle est à effets fixes si le seul effet aléatoire est le terme d’erreur. Il est à effets aléatoires si en dehors de la moyenne générale (qui n’est pas aléatoire) et du terme résiduel, tous les autres termes sont des effets aléatoires. Le modèle est à effets mixtes s’il incorpore à la fois des effets fixes et des effets aléatoires. Le modèle à effets mixtes a donc plusieurs sources de variabilité (effets aléatoires et terme résiduel) et est parfois appelé modèle hiérarchique [113], modèle multi-
niveaux [76] ou encore modèle à composantes de variance [88, 162, 168] (car la variance de $y_{ij}$ se décompose : $\sigma_y^2 = \sigma_\gamma^2 + \sigma_\epsilon^2$).

Dans plusieurs domaines (économie, éducation, sociologie, biologie, etc.), le recours aux données longitudinales a pris une ampleur considérable pour plusieurs raisons dont la principale est la possibilité d'étudier, à la fois, l'aspect dynamique et transversal d'un problème d'intérêt au sein d'une population [64]. En économie par exemple, on utilise l'appellation données de panel pour désigner des données de type longitudinal. Ici, un panel désigne un groupe d'individus suivis de façon répétée au cours du temps. Le terme étude de panel a été utilisé pour la première fois par Lazarsfeld et Fiske [107] qui étudiaient l'effet de la publicité radiodiffusée sur la vente d'un produit commercialisé. Ils entreprirent d'examiner si les personnes achetant leur produit étaient plus susceptibles d'entendre l'annonce. Pour ce faire, ils interrogea à plusieurs reprises un groupe (le panel) de personnes. Considérées comme des séries chronologiques (temporelles) transversales groupées, les données de panel étaient initialement analysées (en économie) par une stratégie en deux étapes :

1) estimer des paramètres transversales en ajustant un modèle (standard) de régression aux données transversales à temps fixé ;

2) utiliser les méthodes propres à l'analyse des séries chronologiques pour modéliser les paramètres transversales estimés.

Bien qu'utile dans certains contextes, cette approche est complètement inappropriées dans d'autres. Par exemple, Edward W. Frees [64] a régressé le taux de divorce (en 1965 et 1975) des cinquante Etats des Etats Unis d'Amérique sur l'aide financière apportées aux familles ayant des enfants à charge. Pour cet exemple, chaque Etat représente une unité statistique observée deux fois (en 1965 et en 1975) au cours du temps, la variable d'intérêt étant le taux de divorce, et l'aide financière, une variable devant aider à mieux comprendre le taux de divorce. Les pentes transversales estimées sont de -0.95% pour 1965 et de -1.0% pour 1975. L'extrapolation de ces pentes négatives donne un résultat complètement différent de celui de l'estimateur de la dynamique du taux de divorce qui est de 2.9%. Afin de contourner ce problème, Theil et Goldberger [184] proposent plutôt une analyse directe (en une seule étape) et fournissent ainsi, en économie, les premières discussions concernant les avantages de l'analyse simultanée de l'aspect transversal et de l'aspect temporel (série chronologique) des données de panel.

Les outils dédiés au traitement de séries chronologiques unidimensionnelles sont bien adaptés à l'analyse de la dynamique des relations entre les observations d'un même sujet. Toutefois, ces méthodes ne sont pas appropriées pour prendre en compte les relations entre différents sujets. Par contre, les méthodes d'analyse dédiées aux séries chronologiques multidimensionnelles peuvent prendre en compte les relations inter-sujets, mais sont utiles dans le cas de données de panel pour un nombre limité de sujets. Aussi, le nombre d'observations par sujet doit-il être au moins de l'ordre de la trentaine afin d'espérer une inférence relativement fiable [64]. Ce qui est vu comme une importante limite du recours aux méthodes de traitement des séries chronologiques pour l'analyse des données de panel. Ainsi, l'approche consistant à considérer les données de panel comme étant des séries chronologiques répétées a très vite montré ses limites et n'a pas connu un grand succès [64].
L'utilisation des données de panel (en économie) est généralement destinée à l'analyse de la dynamique comportementale des sujets au fil du temps. Celle-ci revient à l'analyse de la nature des termes latents de perturbation qui gouvernent les relations (économiques) à l'étude, et qui sont communs aux sujets. Ainsi, le problème central d'intérêt est celui de l'hétérogénéité latente relative aux sujets [133]. L'hétérogénéité entre sujets peut s'expliquer par la dynamique du profil des observations ou par le fait que les observations partagent un même phénomène sous-jacent, non observable, et qui induit une corrélation significative entre les données d'un même sujet. Ainsi, l'hétérogénéité s'interprète ici par le fait que les observations venant du même sujet ont plus tendance à se "ressembler" que celles venant de différents sujets. Finalement, cette hétérogénéité se modélise de deux manières [64]. La première consiste à distinguer les paramètres qui sont communs à tous les sujets de ceux qui sont spécifiques à chaque sujet. Cette approche est connue en économie sous l'appellation de modèle à effets fixes. La seconde méthode ajoute une hypothèse distributionnelle aux paramètres spécifiques aux sujets. Ce qui revient au modèle à effets mixtes ou au modèle à effets aléatoires décrit plus haut.

Bien qu'introduit depuis le 19ᵉ siècle et utilisé dans plusieurs domaines, ce n'est que dans les années 1980 que l'utilité du modèle linéaire à effets mixtes a été fortement mise en lumière par les travaux de Laird et Ware [105] qui en ont proposé une formulation plus générale et plus flexible pouvant naturellement prendre en compte plusieurs particularités des données dont principalement leur aspect déséquilibré(1). La formulation proposée par Laird et Ware [105] est la suivante :

$$y_i = X_i\beta + Z_i\gamma_i + \varepsilon_i, \quad (1.3)$$

avec

$$\gamma_i \sim \mathcal{N}(0, \Gamma), \quad \varepsilon_i \sim \mathcal{N}(0, \Sigma_i) \quad \text{et} \quad \gamma_i \perp \varepsilon_i. \quad (1.4)$$

Dans les équations (1.3) et (1.4), $y_i$ est le vecteur des observations du $i^{\text{e}}$ sujet, $\beta$ appartient à $\mathbb{R}^p$ et est le vecteur des effets fixes, $X_i$ est une matrice de covariables de taille $n_i \times p$, où $n_i$ est le nombre d'observations du $i^{\text{e}}$ sujet. $\gamma_i$ est le vecteur, de taille $q'$, des effets aléatoires relatifs au $i^{\text{e}}$ sujet et $Z_i$ est une matrice de design basée sur des covariables et de taille $n_i \times q'$. $\varepsilon_i$ est le terme résiduel du modèle. $\Gamma$ et $\Sigma_i$ sont des matrices définies positives de taille $q' \times q'$ et $n_i \times n_i$, respectivement. $\Sigma_i$ ne dépend de $i$ que par sa taille. Les sujets étant indépendants les uns des autres, $\gamma_i$ est indépendant de $\gamma_{i'}$ pour tout $i \neq i'$. L'un des avantages de cette formulation est qu'elle ne fait pas de restrictions sur les matrices $X$ et $Z$, et permet une estimation efficace des paramètres du modèle par les méthodes basées sur le maximum de vraisemblance. Ainsi, Laird et Ware [105] ont montré comment l'algorithme EM (Expectation-Maximization) de Dempster, Laird et Rubin [44] peut être utilisé pour estimer les paramètres du modèle dans le cas de données longitudinales. Ce qui représentait, à cette époque, une vraie avancée puisque l'ajustement du modèle, dans un cadre général, à un jeu de données longitudinales posait toujours problème [61]. Peu de temps après, Jennrich et Schluchter [96] proposèrent de recourir à plusieurs méthodes alternatives d'estimation des paramètres du modèle dont principalement l'algorithme de Newton-Raphson et celui des scores de Fisher. Depuis les travaux de Jennrich et Schluchter [96], l'estimation et l'inférence dans le modèle à effets mixtes sont généralement basées sur le maximum de vraisemblance (Maximum

<footnote>(1)Est dit déséquilibré, un jeu de données longitudinales dont les sujets n'ont pas le même nombre d'observations et/ou les observations répétées des sujets n'ont pas été prises à la même date.</footnote>
Likelihood - ML - en anglais) et le maximum de vraisemblance restreinte (Restricted Maximum Likelihood - REML - en anglais). Dans la suite, nous nous focaliserons sur l'estimation des paramètres du modèle linéaire à effets mixtes en dimension *k*, avec $k \ge 1$. Dans la section suivante, nous allons nous consacrer à la revue des méthodes d'estimation des paramètres du modèle en dimension 1.

## 1.2 Modèle linéaire à effets mixtes : différentes méthodes d'estimation des paramètres

Nous allons définir plus explicitement le modèle linéaire à effets mixtes avant de présenter les différentes méthodes d'estimation de ses paramètres.

### 1.2.1 Définition du modèle

Nous faisons ici l'option d'introduire la définition du modèle linéaire à effets mixtes à travers un exemple simple. Supposons que l'on cherche à comprendre l'effet de la qualité nutritionnelle des enfants d'un pays (d'une population) P, pendant leur trois premières années de vie, sur leur croissance.

<table><caption>TABLE 1.1 : Exemple de données longitudinales</caption><thead><tr><th rowspan="2">sujet</th><th rowspan="2">age</th><th rowspan="2">sexe</th><th rowspan="2">MScoreN</th><th colspan="3">Variables d'intérêt</th></tr><tr><th>ScoreN</th><th>poids</th><th>taille</th></tr></thead><tbody><tr><td>1</td><td>0</td><td>F</td><td>63.76</td><td>38.16</td><td>3.14</td><td>47.82</td></tr><tr><td>2</td><td>4</td><td>M</td><td>100.88</td><td>41.46</td><td>4.87</td><td>64.02</td></tr><tr><td>1</td><td>6</td><td>F</td><td>60.98</td><td>41.37</td><td>8.43</td><td>73.21</td></tr><tr><td>3</td><td>0</td><td>M</td><td>93.24</td><td>48.76</td><td>2.82</td><td>44.93</td></tr><tr><td>2</td><td>7</td><td>M</td><td>101.95</td><td>44.79</td><td>8.03</td><td>89.54</td></tr><tr><td>2</td><td>10</td><td>M</td><td>99.24</td><td>48.17</td><td>10.08</td><td>92.14</td></tr><tr><td>1</td><td>16</td><td>F</td><td>85.02</td><td>44.79</td><td>13.96</td><td>86.12</td></tr><tr><td>3</td><td>9</td><td>M</td><td>88.38</td><td>47.91</td><td>8.47</td><td>86.42</td></tr></tbody></table>

On recrute aléatoirement *n* enfants à la naissance que l'on suit jusqu'à l'âge de 36 mois. Tous les mois, on mesure leur poids (en kg), leur taille (en cm), leur score de nutrition (ScoreN) aussi bien que le score de nutrition de leur mère (MScoreN). Pour faire plus simple, intéressons-nous à trois des enfants dont une partie des données est présentée dans le tableau 1.1. L'enfant (le sujet) 1 est une fille qui pèse 3.14 kg et mesure 47.82 cm à la naissance, et a été vue le sixième et le seizième mois, d'après les informations contenues dans le tableau 1.1. Le poids et la taille sont a priori corrélés, et on pourrait commencer par expliquer le poids de l'enfant (le sujet) *i* au *j*ᵉ mois avec le modèle linéaire ci-après :

$$ \text{poids}_{ij} = \beta_0 + \beta_1 \text{taille}_{ij} + \varepsilon_{ij}; \quad i = 1, \dots, n; \quad j = 1, \dots, n_i, \qquad (1.5) $$

où $\beta_0$ est le poids moyen des enfants de la population P, $\beta_1$ mesurant l'effet de la taille sur le poids et $\varepsilon_{ij}$ le terme d'erreur. Les termes poids$_{ij}$ et taille$_{ij}$ désignent respectivement la $j$ᵉ mesure du poids et de la taille du sujet *i*. On pourrait faire l'hypothèse que $\varepsilon_{ij} \sim N(0, \sigma^2)$. Ici, $n = 3$ sujets,
$n_1 = n_2 = 3$ et $n_3 = 2$. Le nombre total d'observations est $N = \sum_{i=1}^{n} n_i = 8$. Puisque les enfants sont nés et grandissent chacun dans des conditions spécifiques, l'idée d'avoir un $\beta_0$ et un effet global $\beta_1$ de la taille sur le poids peut être complètement erronée, et il pourrait être plus pertinent d'introduire dans le modèle, un effet spécifique par sujet comme suit :

$$ \text{poids}_{ij} = \beta_0 + \gamma_i^0 + (\beta_1 + \gamma_i^t)\text{taille}_{ij} + \varepsilon_{ij}; \quad i = 1, \dots, n; \quad j = 1, \dots, n_i. \quad (1.6) $$

Par ailleurs, puisque l'on s'intéresse à tous les enfants de la population P, et non aux seuls n enfants recrutés dans l'échantillon, on fait l'hypothèse :

$$ \begin{pmatrix} \gamma_i^0 \\ \gamma_i^t \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}, \begin{pmatrix} \eta_1^2 & \rho\eta_1\eta_2 \\ \rho\eta_1\eta_2 & \eta_2^2 \end{pmatrix} \right), \quad i = 1, \dots, n. \quad (1.7) $$

Le terme $\gamma_i^0$ (intercepte aléatoire) mesure ainsi l'écart spécifique induit par le sujet i par rapport au poids moyen $\beta_0$ de la population P. $\beta_1$ mesure l'effet moyen global de la taille sur le poids et $\beta_1 + \gamma_i^t$ mesure l'effet spécifique de la taille du sujet i sur son poids ; $\gamma_i^t$ étant la pente aléatoire induit par le sujet i, suivant la variable "taille". En prenant en compte, dans le modèle, le sexe et la qualité nutritionnelle de l'enfant, on a par exemple pour $i = 1, \dots, n$; et $j = 1, \dots, n_i$,

$$ \text{poids}_{ij} = \beta_0 + \gamma_i^0 + (\beta_1 + \gamma_i^t)\text{taille}_{ij} + \gamma_i^s\text{ScoreN}_{ij} + \beta_2\mathbb{1}_{\text{sexe=M}}(\text{sexe}_i) + \beta_3\mathbb{1}_{\text{sexe=F}}(\text{sexe}_i) + \varepsilon_{ij}, \quad (1.8) $$

avec

$$ \gamma_i = \begin{pmatrix} \gamma_i^0 \\ \gamma_i^t \\ \gamma_i^s \end{pmatrix} \sim \mathcal{N}(\mathbf{0}, \Gamma) \quad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2). \quad (1.9) $$

Ce qui s'écrit plus simplement à l'échelle individuelle en regroupant parties fixes et parties aléatoires comme suit :

$$ \text{poids}_i = \beta_0 \mathbb{1}_{n_i} + \beta_1 \text{taille}_i + \beta_2 \mathbb{1}_{\text{sexe=M}}(\text{sexe}_i) \mathbb{1}_{n_i} + \beta_3 \mathbb{1}_{\text{sexe=F}}(\text{sexe}_i) \mathbb{1}_{n_i} \\ + \gamma_i^0 \mathbb{1}_{n_i} + \gamma_i^t \text{taille}_i + \gamma_i^s \text{ScoreN}_i + \varepsilon_i, \quad i = 1, \dots, n, $$

où, $\mathbb{1}_{n_i} = (1, \dots, 1)^\top \in \mathbb{R}^{n_i}$, $\varepsilon_i = (\varepsilon_{i1}, \dots, \varepsilon_{in_i})^\top$, et $\text{taille}_i$ et $\text{ScoreN}_i$ sont des vecteurs qui contiennent les mesures successives du poids et du score de nutrition du sujet i, respectivement. Pour des questions d'identifiabilité du modèle, on met des contraintes sur les paramètres $\beta_2$ et $\beta_3$ qui mesurent l'effet des modalités de la variable "sexe" (variable non numérique) sur le poids, exactement comme dans un modèle linéaire simple par exemple. Si nous faisons l'option d'une
contrainte de type “cellule de référence” sur les modalités du sexe, où la modalité de référence est “sexe = F” (c’est à dire $\beta_3 = 0$), matriciellement, le modèle s’écrit :

$$ \text{poids} = X\beta + Z\gamma + \varepsilon; \quad \gamma \sim \mathcal{N}(\mathbf{0}, \Gamma), \quad \varepsilon \sim \mathcal{N}(\mathbf{0}, \Sigma). \qquad (1.10) $$

Ainsi, $\beta$, $\Gamma$ et $\Sigma$ sont les paramètres (à estimer) du modèle. Ici, $\beta = (\beta_0, \beta_1, \beta_2)^\top$, puisque nous avons fait l'option d'une contrainte d'identifiabilité de type “cellule de référence” sur les modalités du sexe, où la modalité de référence est “sexe = F” dans le modèle exprimé par l'équation (1.10). Sous cette contrainte, $\beta_2$ s'interprète plutôt comme l'effet différentiel du sexe sur le poids, toutes choses étant égales par ailleurs. En considérant le jeu de données fourni par le tableau 1.1, $\gamma = (\gamma_1, \gamma_2, \gamma_3)^\top$ avec $\gamma_i, i = 1, \dots, 3$ définis dans l'équation (1.9). Aussi a-t-on :

$$ X = \begin{pmatrix} 1 & 47.82 & 0 \\ 1 & 64.02 & 1 \\ 1 & 73.21 & 0 \\ 1 & 44.93 & 1 \\ 1 & 89.54 & 1 \\ 1 & 92.14 & 1 \\ 1 & 86.12 & 0 \\ 1 & 86.42 & 1 \end{pmatrix}, $$

$$ Z = \begin{pmatrix} 1 & 47.82 & 38.16 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 64.02 & 41.46 & 0 & 0 & 0 \\ 1 & 73.21 & 41.37 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 44.93 & 48.76 \\ 0 & 0 & 0 & 1 & 89.54 & 44.79 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 92.14 & 48.17 & 0 & 0 & 0 \\ 1 & 86.12 & 44.79 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 86.42 & 47.91 \end{pmatrix}. $$

Si le tableau 1.1 était ordonné par sujet, on aurait pour la matrice $Z$ :

$$ Z = \begin{pmatrix} 1 & 47.82 & 38.16 & . & . & . & . & . & . \\ 1 & 73.21 & 41.37 & . & . & . & . & . & . \\ 1 & 86.12 & 44.79 & . & . & . & . & . & . \\ . & . & . & 1 & 64.02 & 41.46 & . & . & . \\ . & . & . & 1 & 89.54 & 44.79 & . & . & . \\ . & . & . & 1 & 92.14 & 48.17 & . & . & . \\ . & . & . & . & . & . & 1 & 44.93 & 48.76 \\ . & . & . & . & . & . & 1 & 86.42 & 47.91 \end{pmatrix}, $$

où les zéros sont représentés par des points. Le modèle linéaire à effets mixtes peut alors se définir comme suit.

**Définition 1.** Soit $\dagger$ une variable aléatoire numérique, expression d'un attribut descripteur des unités statistiques d'une population, et y un vecteur de $\mathbb{R}^N$ dont les composantes sont des réalisations
de $\dagger$ sur un échantillon de taille $n \le N$. On suppose que $y$ est une réalisation d'un vecteur aléatoire $\mathcal{Y}$. Soit $\gamma$ un vecteur aléatoire de $\mathbb{R}^q$ de réalisation non observable. $X$ et $Z$ deux matrices de dimensions respectives $N \times p$ et $N \times q$. On appelle modèle linéaire à effets mixtes, le modèle statistique défini par :

$$ \gamma \sim \mathcal{P}_{\theta}, \quad \theta \in \Theta, \qquad (1.11) $$

$$ \mathcal{Y}|\gamma \sim \mathcal{Q}_{\theta'}, \quad \theta' \in \Theta', \qquad (1.12) $$

$$ \mathbb{E}[\mathcal{Y}|\gamma] = X\beta + Z\gamma, \qquad (1.13) $$

où $\mathcal{P}_{\theta}$ et $\mathcal{Q}_{\theta'}$ sont des lois de probabilité absolument continues par rapport à la mesure de Lebesgue. $\beta \in \mathbb{R}^p$ est appelé le vecteur des effets fixes et $\gamma$ le vecteur des effets aléatoires.

La quasi-totalité des travaux effectués sur le modèle à effets mixtes (linéaire ou non) choisit la loi gaussienne pour $\mathcal{P}_{\theta}$ [102, 115, 153, 195] pour des raisons de simplicité de calculs, et dans le cas de la version linéaire du modèle, $\mathcal{Q}_{\theta'}$ est aussi choisie comme étant une loi gaussienne [8, 13, 191, 202]. Dans cette thèse, nous nous conformons à ce choix, et notre contribution, relativement à la version multidimensionnelle du modèle, reste dans ce cadre précis. Dans ce cadre, et avec les notations de la définition 1, le modèle s'écrit :

$$ \mathcal{Y} = X\beta + Z\gamma + \varepsilon \qquad (1.14) $$

$$ \gamma \sim \mathcal{N}(\mathbf{0}, \Gamma), \quad \varepsilon \sim \mathcal{N}(\mathbf{0}, \Sigma), \quad \gamma \perp \varepsilon. \qquad (1.15) $$

Le terme résiduel du modèle $\varepsilon$ est supposé indépendant du vecteur des effets aléatoires pour des raisons techniques. $\beta$, $\Gamma$ et $\Sigma$ sont les paramètres à estimer du modèle.

### 1.2.2 Estimation des paramètres du modèle linéaire à effets mixtes

Les méthodes d'estimation des paramètres du modèle peuvent être scindées en deux groupes. Le premier regroupe des "méthodes purement statistiques" et le second les méthodes standards d'optimisation connues en analyse numérique sous l'appellation de méthodes du gradient et incluent : l'algorithme de Newton-Raphson, celui de Marquardt, l'algorithme des scores de Fisher, les méthodes du Quasi-Newton, la descente du gradient, la "coordinate descente" du gradient, etc. Nous nous intéressons ici aux méthodes les plus utilisées du premier groupe que nous exposerons. Les méthodes d'estimation des paramètres du modèle linéaire à effets mixtes que nous exposons ici sont : la méthode dite d'estimation empirique de Bayes, l'algorithme EM (Expectation-Maximization), et celle implémentée par Douglas Bates et ses co-auteurs sous le package lme4 [14] du logiciel R [146]. Cette dernière méthode peut être vue comme une amélioration de la méthode d'estimation empirique de Bayes. La plupart de ces méthodes sont basées sur la maximisation de la vraisemblance des paramètres du modèle. Pour le modèle à effets mixtes, l'estimation des paramètres de variance est restée pendant longtemps un véritable problème. Plusieurs stratégies d'estimation, non basées sur la maximisation de la vraisemblance, proposées dans la littérature ont montré leurs limites [88]. Par exemple, pour un modèle d'analyse de la variance, en présence de données équilibrées, l'estimation des paramètres de variance était faite en posant une égalité
entre les carrés moyens du tableau d'analyse de la variance et leur espérances. Henderson [92] a également développé une méthode analogue, largement utilisée, pour des données déséquilibrées. Progressivement, l'intérêt porté aux méthodes s'appuyant sur la maximisation de la vraisemblance pour l'estimation des paramètres de variance a pris de l'ampleur, car elles présentent plusieurs caractéristiques intéressantes : consistence, normalité asymptotique et efficience (dans le sens décrit par Miller [126]). Par ailleurs, la méthode du maximum de vraisemblance est toujours bien définie, même en présence de contraintes sur certains paramètres de variance. Aussi obtient-on aisément les estimateurs de maximum de vraisemblance et la matrice d'information de Fisher pour n'importe quelle paramétrisation du modèle [88]. Malgré ses intéressantes propriétés, l'approche du maximum de vraisemblance pour l'estimation des composantes de variance est restée longtemps sans être utilisée en pratique puisque, à l'exception des cas simples, le calcul des estimateurs du maximum de vraisemblance exige le recours aux méthodes numériques qui résolvent les problèmes d'optimisation non linéaire à contraintes. Il a fallu attendre l'avènement des ordinateurs électroniques pour surmonter ce problème et accéder à une utilisation plus courante de l'approche. Cependant, même après l'avènement des ordinateurs, il était important d'écrire des algorithmes efficaces pour le calcul des estimateurs de maximum de vraisemblance de composantes de variance.

Dans cette section, nous considérons le modèle linéaire à effets mixtes défini par les équations (1.14) et (1.15), avec $y_{\text{obs}}$ une réalisation observée de $\mathcal{Y}$.

## Estimation empirique de Bayes

La méthode d'estimation empirique de Bayes fait historiquement partie des premières méthodes les plus utilisées pour l'estimation des paramètres du modèle linéaire à effets mixtes. L'approche classique d'estimation est celle du maximum de vraisemblance pour les effets fixes et les paramètres de variance, où les effets aléatoires sont obtenus par l'utilisation de la version étendue du théorème de Gauss-Markov pour l'"estimation" d'effets aléatoires, proposée par Harville [87]. Soit $\theta$ le vecteur de tous les paramètres de variance du modèle. Ces paramètres de variance sont ceux provenant de $\Gamma$ et de $\Sigma$. Les paramètres à estimer sont désormais $\beta$ et $\theta$, en plus d'une réalisation des effets aléatoires conditionnellement à $y_{\text{obs}}$. L'algorithme de l'estimation empirique de Bayes est le suivant :

1. On initialise $\theta$ par une valeur arbitraire $\theta_0$.

2. Avec $\theta$, on écrit la matrice de précision de $y_{\text{obs}}$

$$W = (\Sigma + Z\Gamma Z^{\top})^{-1}. \qquad (1.16)$$

3. L'estimateur de $\beta$ est celui des moindres carrés généralisés

$$\hat{\beta} = (X^{\top} W X)^{-1} X^{\top} W y_{\text{obs}}. \qquad (1.17)$$
4. Une réalisation de $\gamma$, conditionnellement à $y_{\text{obs}}$ s'écrit

$$ \tilde{\gamma} = \mathbb{E} \left[ \gamma | y_{\text{obs}}, \hat{\beta}, \theta \right] = \Gamma Z^{\top} W (y_{\text{obs}} - X \hat{\beta}). \quad (1.18) $$

5. La vraisemblance s'écrit finalement comme fonction uniquement de $\theta$, facile à optimiser. Avec la valeur $\hat{\theta}$ de $\theta$ qui maximise la vraisemblance, on obtient $\hat{\beta}, \hat{\Gamma}$ et $\hat{\Sigma}$.

Cette méthode tire son nom de l'expression de $\tilde{\gamma}$ et n'est en fait qu'une approche unifiée basée sur une combinaison d'estimation par maximum de vraisemblance et de Bayes empirique.

Les estimateurs du maximum de vraisemblance du vecteur des composantes de variance, $\theta$, ne tiennent pas compte de la perte en degré de liberté résultant de l'estimation des effets fixes du modèle, et sont donc biaisés [88, 104]. Patterson et Thompson [139] ont alors proposé l'estimation non biaisée de composantes de variance par la méthode du REML. Cette méthode, vue comme une alternative efficace à l'estimation basée sur la maximisation de la vraisemblance, consiste à maximiser l'intégrale de la vraisemblance ; laquelle intégrale est calculée relativement au vecteur des effets fixes.

## Estimateurs de maximum de vraisemblance au sens de l'algorithme EM

L'algorithme EM [44] est un algorithme statistique itératif qui permet d'approcher plus aisément les estimateurs du maximum de vraisemblance des paramètres d'un modèle statistique qu'on ajuste à un jeu de données incomplètes. Dans le cas du modèle linéaire à effets mixtes, les effets aléatoires sont non observés et les données observées, $y$, peuvent être considérées comme étant des observations incomplètes. Considérons le modèle linéaire à effets mixtes défini par les équations (1.14) et (1.15), et posons $\phi$ le vecteur de tous ses paramètres. À l'itération d'ordre $p+1$, la mise en oeuvre de l'algorithme EM pour calculer la valeur $\phi^{p+1}$ du paramètre à partir de $\phi^p$ se fait en deux étapes comme suit :

**Etape E (Espérance) :** calcul de

$$ Q(\phi|\phi^p) = \mathbb{E}[\log f(y, \gamma)|y, \phi^p], \quad (1.19) $$

où $f$ est la fonction densité de la loi jointe de $y$ et $\gamma$.

**Etape M (Maximisation) :** calcul de $\phi^{p+1}$ tel que

$$ \phi^{p+1} = \arg \max_{\phi} Q(\phi|\phi^p). \quad (1.20) $$

Les estimateurs de maximum de vraisemblance au sens de l'algorithme EM s'écrivent.

$$ \hat{\beta} = (X^{\top}\Sigma^{-1}X)^{-1} X^{\top}\Sigma^{-1} (y - Z\mathbb{E}[\gamma|y, \phi^p]), \quad (1.21) $$

$$ \hat{\Gamma} = \mathbb{V}[\gamma|y, \phi^p] + \mathbb{E}[\gamma|y, \phi^p] \mathbb{E}[\gamma|y, \phi^p]^{\top}, \quad (1.22) $$

et

$$ \hat{\Sigma} = Z\mathbb{V}[\gamma|y, \phi^p]Z^{\top} + (y - X\beta - Z\mathbb{E}[\gamma|y, \phi^p])(y - X\beta - Z\mathbb{E}[\gamma|y, \phi^p])^{\top}, \quad (1.23) $$
où,

$$
\mathbb{E}[\gamma|y, \phi^p] = \mathrm{Cov}(\gamma, y|\phi^p)\mathbb{V}[y|\phi^p]^{-1}(y - \mathbb{E}[y|\phi^p]),
$$

$$
\mathbb{V}[\gamma|y, \phi^p] = \Gamma - \mathrm{Cov}(\gamma, y|\phi^p)\mathbb{V}[y|\phi^p]^{-1}\mathrm{Cov}(\gamma, y|\phi^p)^\top,
$$

$$
\mathbb{E}[y|\phi^p] = X\beta, \quad \mathbb{V}[y|\phi^p] = Z\Gamma Z^\top + \Sigma \text{ et } \mathrm{Cov}(\gamma, y|\phi^p) = \Gamma Z^\top.
$$

Ces estimateurs s'obtiennent simplement en recherchant les points critiques de la fonction

$\phi \mapsto Q(\phi|\phi^p)$, puis en considérant que $\begin{pmatrix} \gamma \\ y \end{pmatrix}$ est un vecteur gaussien.

Estimation des paramètres selon l'approche de Bates et ses co-auteurs [14]

Dans cette section, nous décrivons la stratégie d'optimisation utilisée sous le package lme4 [14] du logiciel R pour obtenir la valeur des estimateurs des paramètres d'un modèle linéaire à effets mixtes, relativement à un jeu de données. Cette stratégie exige que les résidus soient homosédastiques de façon à ce que le modèle s'écrive

$$
(Y|\gamma) \sim \mathcal{N}(X\beta + Z\gamma, \sigma^2 I_n), \qquad (1.24)
$$

$$
\gamma \sim \mathcal{N}(0, \Gamma). \tag{1.25}
$$

Nous rappelons que $X$ est une matrice $n \times p$, $Z$ est une matrice $n \times q$, $I_n$ est la matrice identité d'ordre $n$, $\Gamma$ est une matrice $q \times q$, $\beta \in \mathbb{R}^p$, $\sigma \in \mathbb{R}_+$ et $\gamma \in \mathbb{R}^q$ est le vecteur des effets aléatoires.

Les paramètres du modèle sont désormais : $\beta$, $\sigma$ et $\Gamma$.

$\Gamma$ étant la matrice de variance-covariance des effets aléatoires, elle doit être définie semi-
positive. Il est alors plus pratique d'exprimer le modèle en terme de facteur de covariance relatif,
$\Lambda_{\theta}$ en utilisant la décomposition de Cholesky ci-après

$$
\Gamma = \sigma^2 \Lambda_{\theta} \Lambda_{\theta}^{\top}, \tag{1.26}
$$

où $\theta$ est le vecteur de composantes de variance qui génère $\Lambda_{\theta}$. Bien que $q$, la dimension de $\Gamma$,
soit parfois très grande, la dimension de $\theta$ est généralement faible (moins de 10). L'égalité (1.26)
permet de définir la variable aléatoire $\mathcal{U}$ dénommée *variable sphérique des effets aléatoires* comme
suit

$$
\mathcal{U} \sim \mathcal{N}(0, \sigma^2 I_q), \tag{1.27}
$$

$$
\gamma = \Lambda_{\theta} \mathcal{U}. \tag{1.28}
$$

**Ecriture de la vraisemblance** $y_{\text{obs}}$ étant une réalisation observée de $\mathcal{Y}$, la vraisemblance des paramètres du modèle, conditionnellement à $y_{\text{obs}}$ s'écrit
$$L(\theta, \beta, \sigma^2 | y_{obs}) = \int_{\mathbb{R}^q} f_{y,\mathcal{U}}(y_{obs}, u) du. \qquad (1.29)$$

$f_{y,\mathcal{U}}(y, u)$, la densité jointe de $\mathcal{Y}$ et de $\mathcal{U}$ s'écrit

$$ \begin{align} f_{y,\mathcal{U}}(y, u) &= f_{y|\mathcal{U}}(y|u) f_{\mathcal{U}}(u) \nonumber \\ &= \frac{\exp(-\frac{1}{2\sigma^2} ||y - X\beta - Z\Lambda_{\theta}u||^2)}{(2\pi\sigma^2)^{n/2}} \frac{\exp(-\frac{1}{2\sigma^2} ||u||^2)}{(2\pi\sigma^2)^{q/2}} \nonumber \\ &= \frac{\exp(-\frac{[||y - X\beta - Z\Lambda_{\theta}u||^2 + ||u||^2]}{(2\pi\sigma^2)^{(n+q)/2}}) \tag{1.30} \end{align} $$

Par ailleurs,

$$f_{y,\mathcal{U}}(y_{obs}, u) = f_{\mathcal{U}|\mathcal{Y}}(u|y_{obs}) \int f_{y,\mathcal{U}}(y_{obs}, u) du \qquad (1.31)$$

Ainsi, $f_{y,\mathcal{U}}(y_{obs}, u)$ et $f_{\mathcal{U}|\mathcal{Y}}(u|y_{obs})$ atteignent leur maximum au même point $\mu_{\mathcal{U}|\mathcal{Y}=y_{obs}}$. Et donc,

$$ \begin{align} \left( \begin{array}{c} \mu_{\mathcal{U}|\mathcal{Y}=y_{obs}} \\ \hat{\beta}_{\theta} \end{array} \right) &= \underset{u,\beta}{\arg\max} f_{y,\mathcal{U}}(y_{obs}, u) \nonumber \\ &= \underset{u,\beta}{\arg\min} (||y_{obs} - X\beta - Z\Lambda_{\theta}u||^2 + ||u||^2) \tag{1.32} \end{align} $$

**Remarque** La détermination simultanée de l'espérance conditionnelle des effets aléatoires $\mu_{\mathcal{U}|\mathcal{Y}=y_{obs}}$ (qui est aussi le mode de la distribution conditionnelle des effets aléatoires) et du paramètre $\beta$, prend la forme d'un problème des moindres carrés pénalisés. Cette remarque permet d'obtenir les expressions analytiques de $\mu_{\mathcal{U}|\mathcal{Y}=y_{obs}}$ et de $\beta$, ce qui permet ensuite d'obtenir une écriture expliquite et simple à optimiser de la vraisemblance $L(\theta, \beta, \sigma^2|y_{obs})$ relativement à une valeur fixée de $\theta$. En itérant le processus suivant les valeurs de $\theta$, on obtient l'estimateur $\hat{\theta}$ du maximum de vraisemblance de $\theta$; ce qui fournit l'estimateur de $\Lambda_{\theta}$.

**Notations** Dans la suite, nous adoptons les notations ci-après :

$$g(u, \beta) = \|y_{obs} - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2 \qquad (1.33)$$

$$r_{\theta}^{2} = \min_{u, \beta} (\|y_{obs} - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2) \qquad (1.34)$$

et $\mu_{\mathcal{U}|\mathcal{Y}=y_{obs}} = \mu$ pour raison de simplicité.
$g(u, \beta)$ en fonction de $\mu_{u|Y=y_{obs}}$ et de $\hat{\beta}_{\theta}$ :

$$
\begin{align}
\begin{pmatrix} \mu u | Y=y_{obs} \\ \hat{\beta}_{\theta} \end{pmatrix} &= \operatorname*{arg\,min}_{u,\beta} g(u, \beta) \nonumber \\
&= \operatorname*{arg\,min}_{u,\beta} \left\| \begin{pmatrix} y_{obs} \\ 0 \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix} \begin{pmatrix} u \\ \beta \end{pmatrix} \right\|^2 \tag{1.35}
\end{align}
$$

$$
\Rightarrow \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix}^{\top} \begin{pmatrix} y_{obs} \\ 0 \end{pmatrix} = \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix}^{\top} \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix} \begin{pmatrix} \mu \\ \hat{\beta}_{\theta} \end{pmatrix} \quad \text{Eq. normale} \tag{1.36}
$$

$$
\Leftrightarrow \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix}^{\top} \left[ \begin{pmatrix} y_{obs} \\ 0 \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} & X \\ I_q & 0 \end{pmatrix} \begin{pmatrix} \mu \\ \hat{\beta}_{\theta} \end{pmatrix} \right] = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \quad (1.37)
$$

$$
\Leftrightarrow \begin{cases} (Z\Lambda_{\theta})^{\top} (y_{obs} - Z\Lambda_{\theta}\mu - X\hat{\beta}_{\theta}) + \mu = 0 \\ X^{\top} (y_{obs} - Z\Lambda_{\theta}\mu - X\hat{\beta}_{\theta}) = 0 \end{cases} \tag{1.38}
$$

Par ailleurs,

$$
(1.36) \quad \Rightarrow \quad \begin{pmatrix} (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + I_q & (Z\Lambda_{\theta})^{\top} X \\ X^{\top} Z\Lambda_{\theta} & X^{\top} X \end{pmatrix} \begin{pmatrix} \mu \\ \hat{\beta}_{\theta} \end{pmatrix} = \begin{pmatrix} (Z\Lambda_{\theta})^{\top} y_{obs} \\ X^{\top} y_{obs} \end{pmatrix} \quad (1.39)
$$

On définit le facteur $L_{\theta}$ de Cholesky tel que

$$
L_{\theta} L_{\theta}^{\top} = (Z \Lambda_{\theta})^{\top} Z \Lambda_{\theta} + I_{q}. \qquad (1.40)
$$

$g$ est de classe $C^{\infty}$. En posant $x = \begin{pmatrix} u \\ \beta \end{pmatrix}$ et $\hat{x} = \begin{pmatrix} \mu \\ \hat{\beta}_{\theta} \end{pmatrix}$, on a :

$$
g(x) = g(\hat{x}) + g'(\hat{x})(x - \hat{x}) + \frac{1}{2}(x - \hat{x})^{\top} g''(\hat{x})(x - \hat{x}) + \dots + \dots \quad (1.41)
$$

$$
\frac{\partial g}{\partial u}(\hat{x}) = 2L_{\theta}L_{\theta}^{\top}(\mu - \mu) = 0; \quad \frac{\partial g}{\partial \beta}(\hat{x}) = -2X^{\top}(y_{obs} - Z\Lambda_{\theta}\mu - X\hat{\beta}_{\theta}) = 0 \text{ d'après (1.38)}.
$$

Ainsi, $g'(\hat{x}) = (0,0)$.

$$
\frac{\partial g}{\partial u^2}(x) = 2L_{\theta}L_{\theta}^{\top}; \quad \frac{\partial g}{\partial u \partial \beta}(x) = 2X^{\top}Z\Lambda_{\theta}; \quad \frac{\partial g}{\partial \beta \partial u}(x) = 2(Z\Lambda_{\theta})^{\top}X; \quad \frac{\partial g}{\partial \beta^2}(x) = 2X^{\top}X
$$

Ainsi,

$$
g''(\hat{x}) = 2 \begin{pmatrix} L_{\theta} L_{\theta}^{\top} & (Z \Lambda_{\theta})^{\top} X \\ Z \Lambda_{\theta} X^{\top} & X^{\top} X \end{pmatrix} = 2M \quad (1.42)
$$

$$
\text{et } g^{(n)}(x) = 0 \quad \forall n \ge 3.
$$

D'où
$$
\begin{align}
g(x) &= g(\hat{x}) + (x - \hat{x})^{\top} M (x - \hat{x}) \nonumber \\
     &= r_{\theta}^{2} + (x - \hat{x})^{\top} M (x - \hat{x}) \tag{1.43}
\end{align}
$$

En faisant la décomposition de Cholesky de $M$ ci-après

$$
M = \begin{pmatrix} L_{\theta}^{\top} & R_{ZX} \\ 0 & R_{X} \end{pmatrix}^{\top} \begin{pmatrix} L_{\theta}^{\top} & R_{ZX} \\ 0 & R_{X} \end{pmatrix}, \quad (1.44)
$$

on a :

$$
\begin{align}
g(u, \beta) &= r_{\theta}^{2} + \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix}^{\top} \begin{pmatrix} L_{\theta} & 0 \\ R_{ZX}^{\top} & R_{X}^{\top} \end{pmatrix} \begin{pmatrix} L_{\theta}^{\top} & R_{ZX} \\ 0 & R_{X} \end{pmatrix} \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix} \nonumber \\
&= r_{\theta}^{2} + \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix}^{\top} \left[ \begin{pmatrix} L_{\theta} \\ R_{ZX}^{\top} \end{pmatrix} \begin{pmatrix} L_{\theta}^{\top} & R_{ZX} \\ 0 & R_{X} \end{pmatrix} + \begin{pmatrix} 0 \\ R_{X}^{\top} \end{pmatrix} \begin{pmatrix} 0 & R_{X} \end{pmatrix} \right] \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix} \nonumber \\
&= r_{\theta}^{2} + \left\| \begin{pmatrix} L_{\theta}^{\top} & R_{ZX} \end{pmatrix} \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix} \right\|^{2} + \left\| \begin{pmatrix} 0 & R_{X} \end{pmatrix} \begin{pmatrix} u - \mu \\ \beta - \hat{\beta}_{\theta} \end{pmatrix} \right\|^{2} \tag{1.45} \\
&= r_{\theta}^{2} + \left\| L_{\theta}^{\top} (u - \mu) + R_{ZX} (\beta - \hat{\beta}_{\theta}) \right\|^{2} + \left\| R_{X} (\beta - \hat{\beta}_{\theta}) \right\|^{2} \tag{1.46}
\end{align}
$$

Expression de la vraisemblance à $\theta$ fixé

$$
\begin{align}
L(\theta, \beta, \sigma^2 | y_{obs}) &= \int_{\mathbb{R}^q} \frac{\exp(g(u, \beta))}{(2\pi\sigma^2)^{(n+q)/2}} du \nonumber \\
&= \frac{\exp\left[-\frac{r_{\theta}^2 + \|R_X(\beta - \hat{\beta}_{\theta})\|^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{n/2}} \int_{\mathbb{R}^q} \frac{1}{(2\pi\sigma^2)^{q/2}} \exp\left[-\frac{\|L_{\theta}^{\top}(u - \mu) + R_{ZX}(\beta - \hat{\beta}_{\theta})\|^2}{2\sigma^2}\right] du \tag{1.47}
\end{align}
$$

En posant $v = L_{\theta}^{\top}(u - \mu) + R_{ZX}(\beta - \hat{\beta}_{\theta})$, on a :

$$
\begin{align*}
dv &= \left| \frac{Dv}{Du} \right| du \\
   &= \left| L_{\theta}^{\top} \right| du
\end{align*}
$$

et
$$
\begin{align}
L(\theta, \beta, \sigma^2 | y_{obs}) &= \frac{\exp\left[-\frac{r_{\theta}^2 + \|R_X(\beta - \hat{\beta}_{\theta})\|^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{n/2}|L_{\theta}|} \underbrace{\int_{\mathbb{R}^q} \frac{1}{(2\pi\sigma^2)^{q/2}} \exp\left[-\frac{\|v\|^2}{2\sigma^2}\right] dv}_{=1} \nonumber \\
&= \frac{\exp\left[-\frac{r_{\theta}^2 + \|R_X(\beta - \hat{\beta}_{\theta})\|^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{n/2}|L_{\theta}|}. \tag{1.48}
\end{align}
$$

**Estimateur de la déviance du modèle à θ fixé** $\hat{\beta}_{\theta}$ est l'estimateur du maximum de vraisem-
blance de $\beta$ à $\theta$ fixé, et en remplaçant $\beta$ par $\hat{\beta}_{\theta}$ dans l'égalité (1.48), on obtient l'estimateur
ci-après de la vraisemblance

$$
\tilde{L}(\theta, \sigma^2 | y_{obs}) = \frac{\exp\left[-\frac{r_{\theta}^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{n/2} |L_{\theta}|}. \quad (1.49)
$$

Ce qui fournit l'estimateur de la déviance ci-après

$$
\begin{align}
\tilde{l}(\theta, \sigma^2 | y_{obs}) &= -2 \log \left( \tilde{L}(\theta, \sigma^2 | y_{obs}) \right) \nonumber \\
&= \frac{r_{\theta}^2}{\sigma^2} + n \log(2\pi\sigma^2) + \log(|L_{\theta}|^2). \tag{1.50}
\end{align}
$$

Vu comme fonction de $\sigma^2$, $\tilde{l}(\theta, \sigma^2|y_{obs})$ est continument dérivable et on a :

$$
\begin{gather*}
\frac{\partial(\tilde{l}(\theta, \sigma^2 | y_{obs}))}{\sigma^2} = -\frac{r_{\theta}^2}{\sigma^4} + \frac{n}{\sigma^2} \\
\frac{\partial(\tilde{l}(\theta, \sigma^2 | y_{obs}))}{\sigma^2} = 0 \iff \sigma^2 = \frac{r_{\theta}^2}{n}.
\end{gather*}
$$

Ainsi,

$$
\hat{\sigma}^2 = \frac{r_{\theta}^2}{n} \tag{1.51}
$$

est l'estimateur du maximum de vraisemblance de $\sigma^2$ à $\theta$ fixé. Ce qui fournit, en remplaçant $\sigma^2$ par $\hat{\sigma}^2$ dans l'égalité (1.50), l'estimateur de la déviance du modèle à $\theta$ fixé ci-après :

$$
\tilde{l}(\theta|y_{obs}) = \log(|L_{\theta}|^2) + n \left[ 1 + \log \left( \frac{2\pi r_{\theta}^2}{n} \right) \right]. \quad (1.52)
$$

Finalement,

$$
\hat{\theta} = \arg \min_{\theta} \tilde{l}(\theta|y_{obs}) \quad (1.53)
$$
**Récapitulatif de la méthode** La méthode d'estimation des paramètres d'un modèle linéaire mixte, utilisée sous le package lme4 du logiciel R se résume en ces quelques étapes qui suivent :

1. On fixe une valeur du paramètre de composante de variance $\theta$, puis on engendre $\Lambda_{\theta}$

2. on calcule $L_{\theta}$ solution de l'équation (1.40)

3. on calcule $\hat{\beta}_{\theta}$ et $\mu$ solutions de l'équation (1.39)

4. on calcule $r_{\theta}^2 = \|y_{obs} - X\hat{\beta}_{\theta} - Z\Lambda_{\theta}\mu\|^2 + \|\mu\|^2$ issu des équations (1.32) et (1.34)

5. on calcule $\tilde{l}(\theta|y_{obs})$ dont l'expression est donnée par l'égalité (1.52)

6. on reprend les opérations 1. à 5. jusqu'à retrouver la valeur $\hat{\theta}$ de $\theta$ qui minimise $\tilde{l}(\theta|y_{obs})$.

$\hat{\beta}_{\theta}$ est l'estimateur du maximum de vraisemblance des effets fixes, et $\tilde{\gamma} = \Lambda_{\theta}\mu$ est le mode de la distribution des effets aléatoires, conditionnellement à $y_{obs}$. $\tilde{\gamma}$ est considéré comme étant la réalisation de $\gamma$, conditionnellement à $y_{obs}$. $\hat{\sigma}^2$ est l'estimateur du maximum de vraisemblance du paramètre de variance résiduelle, tandis que $\hat{\Gamma} = \hat{\sigma}^2\Lambda_{\theta}\Lambda_{\theta}^{\top}$ est l'estimateur du maximum de vraisemblance de la matrice de variance-covariance des effets aléatoires.

**Critère REML pour l'estimation non biaisée des composantes de variance** Le critère REML (Restricted Maximum Likelihood) est la vraisemblance basée sur les composantes de variance et permet d'obtenir un estimateur non biaisé de ces composantes. Il s'écrit

$$L_R(\theta, \sigma^2|y_{obs}) = \int_{\mathbb{R}^p} L(\theta, \beta, \sigma^2|y_{obs})d\beta. \quad (1.54)$$

En utilisant l'égalité (1.48) et en posant $t = R_X(\beta - \hat{\beta}_{\theta})$, on a $dt = |R_X|d\beta$ et

$$L_R(\theta, \sigma^2|y_{obs}) = \frac{1}{|L_{\theta}|} \int_{\mathbb{R}^p} \frac{\exp\left[-\frac{r_{\theta}^2 + \|R_X(\beta - \hat{\beta}_{\theta})\|^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{n/2}} d\beta \quad (1.55)$$

$$= \frac{\exp(-r_{\theta}^2/(2\sigma^2))}{|L_{\theta}| |R_X| (2\pi\sigma^2)^{(n-p)/2}} \underbrace{\int_{\mathbb{R}^p} \frac{\exp\left[-\frac{\|t\|^2}{2\sigma^2}\right]}{(2\pi\sigma^2)^{p/2}} dt}_{=1} \quad (1.56)$$

$$\begin{align*} \frac{\partial (-2 \log(L_R))}{\partial (\sigma^2)} &= 0 &\iff & -\frac{r_{\theta}^2}{\sigma^4} + (n-p)\frac{1}{\sigma^2} &= 0 \\ &\implies & \hat{\sigma}^2 R &= \frac{r_{\theta}^2}{n-p} \tag{1.57} \end{align*}$$

Et donc,
$$
\begin{align}
\tilde{l}_R(\theta|y_{obs}) &= -2(\log(L_R)) \nonumber \\
&= \log(|L_\theta|^2) + \log(|R_X|^2) + (n-p) \left[ 1 + \log\left(\frac{2\pi r_\theta^2}{n-p}\right) \right] \tag{1.58}
\end{align}
$$

puis,

$$
\hat{\theta}_R = \arg \min_{\theta} \tilde{l}_R(\theta|y_{obs}). \tag{1.59}
$$

1.3 Contributions méthodologiques de la thèse

Dans cette thèse, nos contributions sont de deux ordres. Le premier est relatif à l'estimation des
paramètres de la version multidimensionnelle du modèle linéaire à effets mixtes, et comprend deux
articles scientifiques qui représentent les chapitres 2 et 3 du présent manuscrit. Le second quant
à lui s'appuie sur une version modifiée de la méthode d'estimation utilisée dans le package lme4
pour proposer une procédure de sélection d'effets fixes dans le modèle linéaire à effets mixtes.

En ce qui concerne la définition de la version multidimensionnelle du modèle, nous pensons,
pour des raisons de simplicité et pour faciliter la lecture, qu’il serait plus utile de se limiter à
la définition de la version bi-dimensionnelle, puisqu’une généralisation en dimension supérieur
devient naturelle et intuitive. En dimension $k = 2$, le modèle s’écrit

$$
\begin{align}
y_1 &= X_1\beta_1 + Z_1\gamma_1 + \varepsilon_1 \notag \\
y_2 &= X_2\beta_2 + Z_2\gamma_2 + \varepsilon_2 \tag{1.60}
\end{align}
$$

$$
\gamma = \begin{pmatrix} \gamma_1 \\ \gamma_2 \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \mathbf{\Gamma} = \begin{pmatrix} \Gamma_1 & \Gamma_{12} \\ \Gamma_{12}^\top & \Gamma_2 \end{pmatrix} \right), \quad (1.61)
$$

$$
\epsilon = \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \mathbf{\Sigma} = \begin{pmatrix} \Sigma_1 & \Sigma_{12} \\ \Sigma_{12}^\top & \Sigma_2 \end{pmatrix} \right); \quad \gamma \perp \epsilon. \qquad (1.62)
$$

Pour $k \in \{1, 2\}$, $\beta_k$ et $\gamma_k$ désignent respectivement le vecteur des effets fixes et le vecteur des effets aléatoires, et $\varepsilon_k$ est la composante résiduelle. $X_k$ est une matrice de covariables et $Z_k$ une matrice de design faite de covariables. $\dim(X_k) = N \times p_k$ et $\dim(Z_k) = N \times q_k$, avec N le nombre total d'observations composant le jeu de données auquel le modèle est ajusté. $p_k$ et $q_k$ sont respectivement le nombre d'effets fixes et le nombre d'effets aléatoires dans la dimension k du modèle. Les symboles en gras représentent les paramètres qui ont des composantes dans plus d'une dimension du modèle. Par exemple, $\Sigma_1$ intervient seulement dans la dimension 1 du modèle alors que $\Sigma$ intervient dans les deux dimensions du modèle de par ses composantes.

Dans cette thèse, nos contributions, relativement à la version multidimensionnelle du modèle
ont été faites sous l'hypothèse d'indépendance entre les variables réponses conditionnellement
aux effets aléatoires. C'est à dire que les résidus de deux différentes dimensions du modèle sont indépendants. Ce qui se traduit dans la version bi-dimensionnelle par

$$ \Sigma_{12} = 0. \qquad (1.63) $$

Généralement dans la littérature, cette hypothèse est faite pour des raisons techniques relatives au calcul des estimateurs des paramètres de variance résiduelle (voir par exemple, [172], [161], [56] et [6]). Aussi peut-on toujours faire cette hypothèse parce qu'on aura choisi de concentrer toute la dépendance entre réponses marginales dans les effets aléatoires.

Dans nos travaux, nous avons procédé à une revue de littérature des différentes stratégies d'analyse et modèles ajustés aux données longitudinales multidimensionnelles. En effet, les différentes approches utilisées dans la littérature peuvent être scindées en deux groupes à savoir : 1) celles qui spécifient la distribution jointe de toutes les $m$ variables de réponse $y_1, \dots, y_m$ (également noté sous la forme d'un vecteur $\mathbf{y} = (y_1^\top, \dots, y_m^\top)^\top$) sans recourir à des structures latentes et 2) celles qui font une modélisation basée sur des structures latentes.

Dans le premier groupe, on retrouve trois différentes approches. D'abord une spécification directe de la structure de corrélation de $\mathbf{y}$ où l'on peut citer par exemple les travaux de Galecki [68] et ceux de O'Brien [135] qui ont factorisé la matrice de covariance de $\mathbf{y}$ en se servant du produit de Kronecker afin d'avoir des modèles plus parcimonieux dans le cas d'analyse de données longitudinales équilibrées. D'autres modèles comme les copules [132, 175] utilisés par exemple dans les travaux de Lambert [106] peuvent être aussi classés dans ce premier sous-groupe d'approches au même titre que certains modèles de séries chronologiques comme les ARIMA (Auto Regressive Integrated Moving Average) et les VAR (Vector Auto Regressive) utilisés par exemple par MaCurdy [121] et Tschacher [188] dans leurs travaux. Ensuite une analyse sans la modélisation explicite de la structure de corrélation de $\mathbf{y}$ où l'on peut principalement citer les GEE (Generalized Estimating Equations) proposés par Liang [110] en tant qu'extension du modèle linéaire généralisé pour l'analyse de données longitudinales. Des exemples d'utilisation des GEE dans le contexte d'analyse de données longitudinales multidimensionnelles incluent les travaux de Prennice [143], ceux de Liang [111] et de Rochon [156]. Enfin dans ce premier groupe, retrouve-t-on aussi les modèles d'analyses marginales conditionnelles où l'on évite une spécification directe de la distribution jointe de $\mathbf{y}$ (que nous notons $f(\mathbf{y})$) en la factorisant (par exemple en dimension 2, on peut bien écrire $f(y_1, y_2) = f(y_1|y_2)f(y_2)$ où $y_2$ joue le rôle de covariable).

Le second groupe d'approches quant à lui se subdivise en deux sous-groupes dont le premier consiste en une réduction de la dimension de $\mathbf{y}$ et le second, aux modèles connus sous le nom de modèles à effets mixtes. Les stratégies basées sur la réduction de la dimension de $\mathbf{y}$ consistent à utiliser une méthode de type analyse en composantes principales afin de résumer l'information contenue dans $\mathbf{y}$ en un petit nombre de facteurs latents pertinents auxquels il est ensuite ajusté des modèles unidimensionnels standards et indépendants (voir par exemple [49, 123, 137]). En ce qui concerne le modèle à effets mixtes, plusieurs auteurs dont Reinsel [151], Maccallum [120] et Beckett [18] ont proposé son utilisation dans le cadre d'analyse de données longitudinales multidimensionnelles. Bien qu'apparaissant comme étant le modèle le plus souvent utilisé en matière d'analyse de données longitudinales et plus généralement de données multi-niveaux, son utilisation dans le cas multidimensionnel se heurte à des problèmes numériques souvent dus à une dimension
trop élevée du vecteur de tous les effets aléatoires marginaux. Pour contourner ce problème, Fieuws et ses co-auteurs [57] ont proposé comme estimateurs des paramètres du modèle joint, la moyenne des estimateurs issus de tous les $m(m-1)/2$ modèles bi-dimensionnels. Toutes ces approches, leurs forces aussi bien que leurs faiblesses ont été amplement discutées avec davantage de détails dans le deuxième chapitre du présent manuscrit.

Les outils logiciels dédiés à l'ajustement de modèles linéaires à effets mixtes multidimensionnels sont vraiment rares ; ce qui fait que parfois, certains auteurs organisent les données de façon à pouvoir utiliser des outils initialement destinés aux analyses unidimensionnelles. C'est le cas par exemple des travaux de Thiébaut et ses co-auteurs [185] qui ont utilisé la procédure Proc MIXED du logiciel SAS afin d'ajuster un modèle bi-dimensionnel à effets mixtes à des profils longitudinaux de biomarqueurs relatifs au VIH SIDA. Récemment, le programme SabreR [39] du logiciel R, destiné à ajuster des modèles linéaires mixtes à au plus trois dimensions et avec seulement des intercepts aléatoires a été enlevé du dépôt. Schafer et Yucel [161] dont nous avons également mentionné les travaux au chapitre 2, ont mis au point un programme de calculs nommé mlmmm [203] basé sur leur méthode, sous le logiciel R.

Dans cette thèse, nos contributions sont relatives à l'utilisation du modèle linéaire à effets mixtes multidimensionnel et à une amélioration de la qualité de l'estimation de ses paramètres.

Après une bonne lecture de la littérature, on peut remarquer que l'estimation des paramètres du modèle linéaire à effets mixtes multidimensionnel se fait le plus souvent par l'algorithme EM, surtout dans le cadre d'analyse de données longitudinales. C'est pour cela que dans un premier temps, en utilisant l'algorithme EM, nous proposons des estimateurs écrits sous une forme plus générale que celle proposée dans la littérature (voir par exemple [6], [161] et [172]). Ces estimateurs ont une forme plus générale en ce sens qu'ils sont utilisables à la fois pour des données longitudinales et des données multi-niveaux. Par la suite, nous proposons une version multidimensionnelle de la méthode d'estimation proposée par Bates et ses co-auteurs [16]. Bien que restreinte au cas où les résidus dimensionnels sont homoscédastiques et indépendants des effets aléatoires, cette stratégie d'estimation s'est montrée plus robuste par rapport au point de départ (de l'algorithme) et fournit des estimations plus consistantes (sur nos simulations) que celles obtenues avec l'algorithme EM, surtout en ce qui concerne les paramètres de variance des effets aléatoires.

Par ailleurs, nous proposons une procédure de sélection d'effets fixes dans la version unidimensionnelle du modèle. Laquelle procédure est basée sur une pénalisation itérée de type $L_2$ de la log-vraisemblance profilée de certains paramètres du modèle en vue d'approcher les performances d'une pénalité $L_0$. Il importe ici de signaler que ce travail est le prélude d'un autre à venir et qui consistera à faire de la sélection de variables dans la version multidimensionnelle du même modèle. De façon un peu plus détaillée, nous présentons dans les paragraphes qui suivent un résumé de nos différentes contributions.

## Contributions relatives au chapitre 2 : Analyse longitudinale multi-variée et test de corrélation bi-variée

Dans le contexte de l'estimation des paramètres du modèle linéaire multidimensionnel à effets mixtes, nous proposons des expressions plus générales des estimateurs en recourant à l'algorithme EM, exactement comme nous avons commencé par le faire en dimension $k = 1$ (voir
Section 1.2.2). Ces estimateurs ont des expressions plus générales en ce sens qu'ils sont aussi uti-
lisables tels quels dans le cadre plus général de l'estimation des paramètres d'un modèle linéaire
multidimensionnel multi-niveaux. Aussi est-il important de rappeler qu'un jeu de données longi-
tudinales peut-être vu comme étant un jeu de données multi-niveaux. Ces estimateurs du maxi-
mum de vraisemblance obtenus via l'algorithme EM s'écrivent en dimension 2 comme suit. Pour
$k \in \{1, 2\}$,

$$
\hat{\beta}_k = (X_k^\top \Sigma_k^{-1} X_k)^{-1} X_k^\top \Sigma_k^{-1} (y_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \boldsymbol{\phi}^p]), \quad (1.64)
$$

$$
\hat{\Gamma} = \mathbb{V}[\gamma | \mathbf{y}, \phi^p] + \mathbb{E}[\gamma | \mathbf{y}, \phi^p] \mathbb{E}[\gamma | \mathbf{y}, \phi^p]^\top, \quad (1.65)
$$

$$
\hat{\Sigma}_k = Z_k \mathbb{V}[\gamma_k | \mathbf{y}, \phi^p] Z_k^\top + (y_k - X_k \beta_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \phi^p])(y_k - X_k \beta_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \phi^p])^\top, \quad (1.66)
$$

où,

$$
\mathbb{E}[\gamma_k | \mathbf{y}, \phi^p] = \mathrm{Cov}(\gamma_k, \mathbf{y} | \phi^p) \mathbb{V}[\mathbf{y} | \phi^p]^{-1} (\mathbf{y} - \mathbb{E}[\mathbf{y} | \phi^p]), \quad (1.67)
$$

$$
\mathbb{V}[\gamma_k | \mathbf{y}, \phi^p] = \mathbf{\Gamma}_k - \mathrm{Cov}(\gamma_k, \mathbf{y}|\phi^p) \mathbb{V}[\mathbf{y}|\phi^p]^{-1} \mathrm{Cov}(\gamma_k, \mathbf{y}|\phi^p)^{\top} \quad (1.68)
$$

et

$$
\mathbb{V}[\mathbf{y}|\phi^p] = \begin{pmatrix} Z_1 \Gamma_1 Z_1^\top + \Sigma_1 & Z_1 \Gamma_{12} Z_2^\top \\ Z_2 \Gamma_{12}^\top Z_1^\top & Z_2 \Gamma_2 Z_2^\top + \Sigma_2 \end{pmatrix}, \quad \mathrm{Cov}(\gamma, \mathbf{y}|\phi^p) = \begin{pmatrix} Z_1 \Gamma_1 & Z_1 \Gamma_{12} \\ Z_2 \Gamma_{12}^\top & Z_2 \Gamma_2 \end{pmatrix}^\top, \tag{1.69}
$$

$$
\mathrm{Cov}(\gamma_1, \mathbf{y}|\phi^p) = \begin{pmatrix} Z_1 \Gamma_1 \\ Z_2 \Gamma_{12}^\top \end{pmatrix}^\top, \quad \mathrm{Cov}(\gamma_2, \mathbf{y}|\phi^p) = \begin{pmatrix} Z_1 \Gamma_{12} \\ Z_2 \Gamma_2 \end{pmatrix}^\top. \qquad (1.70)
$$

Nous avons procédé à une investigation empirique de la précision de ces estimateurs en calcu-
lant l’erreur quadratique moyenne de plusieurs estimations faites sur 1000 replications de jeux de
données longitudinales simulées sous le modèle bi-dimensionnel ci-contre. En identifiant les sujets
par *i*, et pour *i* = 1, ..., *n*

$$
\begin{align}
y_{1i} &= X_{1i}\beta_1 + Z_{1i}\gamma_{1i} + \varepsilon_{1i} \nonumber \\
y_{2i} &= X_{2i}\beta_2 + Z_{2i}\gamma_{2i} + \varepsilon_{2i}, \tag{1.71}
\end{align}
$$

avec

$$
\gamma_i = \begin{pmatrix} \gamma_{1i} \\ \gamma_{2i} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \bar{\mathbf{\Gamma}} = \begin{pmatrix} \bar{\Gamma}_1 & \bar{\Gamma}_{12} \\ \bar{\Gamma}_{12}^\top & \bar{\Gamma}_2 \end{pmatrix} \right) \quad (1.72)
$$
<table><caption>TABLE 1.2 : Erreur quadratique moyenne des EM-estimateurs avec l'intervalle de confiance à 95% calculés sur 1000 replications pour différentes valeurs du nombre *n* de sujets et du nombre *N* d'observations.</caption><thead><tr><th>Paramètre</th><th>n</th><th>N = 600</th><th>N = 1000</th><th>N = 3000</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50</td><td>1.27 (0.01 - 4.64)</td><td>1.10 (0.00 - 3.83)</td><td>0.71 (0.00 - 2.57)</td></tr><tr><td>60</td><td>1.38 (0.00 - 5.56)</td><td>0.99 (0.00 - 3.90)</td><td>0.65 (0.00 - 2.44)</td></tr><tr><td>100</td><td>1.31 (0.00 - 4.79)</td><td>0.87 (0.00 - 3.41)</td><td>0.46 (0.00 - 1.58)</td></tr><tr><td>300</td><td>1.64 (0.00 - 5.93)</td><td>0.80 (0.00 - 3.23)</td><td>0.29 (0.00 - 1.13)</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>50</td><td>2.13 (0.01 - 7.22)</td><td>1.73 (0.00 - 6.21)</td><td>1.08 (0.00 - 4.17)</td></tr><tr><td>60</td><td>2.12 (0.00 - 8.06)</td><td>1.63 (0.01 - 6.66)</td><td>0.95 (0.00 - 3.52)</td></tr><tr><td>100</td><td>1.88 (0.00 - 7.14)</td><td>1.29 (0.00 - 4.76)</td><td>0.62 (0.00 - 2.24)</td></tr><tr><td>300</td><td>2.30 (0.00 - 8.54)</td><td>1.23 (0.00 - 4.57)</td><td>0.43 (0.00 - 1.65)</td></tr><tr><td rowspan="4">σ<sub>1</sub></td><td>50</td><td>0.03 (0.00 - 0.10)</td><td>0.10 (0.00 - 0.07)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.07)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>100</td><td>0.04 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.08)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>300</td><td>0.06 (0.00 - 0.22)</td><td>0.03 (0.00 - 0.11)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td rowspan="4">σ<sub>2</sub></td><td>50</td><td>0.05 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>60</td><td>0.06 (0.00 - 0.21)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>100</td><td>0.06 (0.00 - 0.24)</td><td>0.04 (0.00 - 0.15)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>300</td><td>0.09 (0.00 - 0.36)</td><td>0.04 (0.00 - 0.18)</td><td>0.01 (0.00 - 0.05)</td></tr><tr><td rowspan="4">Γ̄</td><td>50</td><td>400.91 (1.52 - 1274.32)</td><td>670.46 (2.36 - 2536.73)</td><td>489.07 (2.02 - 1840.60)</td></tr><tr><td>60</td><td>706.45 (2.58 - 2497.42)</td><td>620.34 (2.70 - 2293.58)</td><td>408.79 (0.85 - 1597.77)</td></tr><tr><td>100</td><td>701.50 (3.79 - 2747.70)</td><td>477.55 (2.08 - 1839.26)</td><td>283.25 (1.97 - 1023.77)</td></tr><tr><td>300</td><td>798.28 (3.65 - 2603.65)</td><td>547.83 (3.03 - 1721.46)</td><td>199.54 (0.81 - 736.08)</td></tr></tbody></table>

et

$$ \varepsilon_i = \begin{pmatrix} \varepsilon_{1i} \\ \varepsilon_{2i} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \bar{\Sigma} = \begin{pmatrix} \sigma_1^2 \mathbf{I}_{N_{1i}} & 0 \\ 0 & \sigma_2^2 \mathbf{I}_{N_{2i}} \end{pmatrix} \right) \quad (1.73) $$

$N_{1i}$ et $N_{2i}$ sont les dimensions de $y_{1i}$ et $y_{2i}$, respectivement.

Les résultats des estimations faites sur ces jeux de données de différentes tailles ($n \in \{50, 60, 100, 300\}$, $N \in \{600, 1000, 3000\}$) prouvent empiriquement que les estimateurs des effets fixes et des paramètres de variance résiduelle s'estiment beaucoup mieux (asymptotiquement consistants) que les paramètres de variance des effets aléatoires (voir Tableau 1.2). Ces résultats suggèrent aussi, en ce qui concerne l'analyse de données longitudinales, un bon compromis entre le nombre *n* de sujets et le nombre *N* d'observations en vue d'avoir des estimations d'une qualité relativement acceptable, surtout lorsque *N* n'est pas très grand.

Dans ce chapitre qui a fait l'objet d'une publication scientifique dans PloS one (voir [1]), nous avons également construit un test de rapport de vraisemblance basé sur ces EM-estimateurs afin de tester la significativité de la dépendance entre deux dimensions du modèle ($y_k$ et $y_{k'}$, $k \neq k'$). Dans le cadre de l'analyse de données longitudinales, ce test de corrélation bi-varié a montré de bonnes performances empiriques avec des données simulées. Ce test peut en effet détecter une corrélation de signal faible ($\rho = 0.3$ avec un AUC= 0.81, voir Figure 1.1), même avec des jeux de données de taille modeste ($n = 60$ sujets, avec $N = 600$ observations). En utilisant ce test de corrélation bi-varié, nous avons illustré, par des simulations, certaines des conséquences du choix d'analyser séparément plutôt que conjointement deux variables réponses. Deux variables dépendantes qui sont prouvées non significativement corrélées après le test, seront donc analysées
<figure><img src="image_2.png" /><figcaption>FIGURE 1.1 : Analyse empirique de la puissance du test de corrélation bi-varié. Les valeurs des AUC des courbes ROC avec leurs intervalles de confiance sont calculées, utilisant différentes valeurs du coefficient de corrélation linéaire ρ, du nombre de sujets (n) et du nombre d'observations (N). Le panel de gauche pour N = 600 et n ∈ {50, 60, 100}. Le panel de droite pour N = 3000 et n ∈ {50, 100, 300}.</figcaption></figure>

séparément. Nous pensons que cette stratégie pourrait être utilisée dans une procédure de sélection pas-à-pas du modèle conjoint faisant intervenir toutes les variables dépendantes dans le but d'avoir un modèle conjoint plus parcimonieux en terme de nombre de paramètres de variance à estimer pour les effets aléatoires.

Contributions relatives au chapitre 3 : Estimateurs consistants pour les paramètres du modèle linéaire multidimensionnel à effets mixtes

Dans le chapitre 2, nous avons remarqué que les EM-estimateurs des composantes de variance des effets aléatoires sont ceux qui ont des erreurs quadratiques les plus élevées, même dans le cas d'homoscédasticité des résidus marginaux (dimensionnels) où nous sommes restés pour simuler et estimer. Dans le chapitre 3, nous restons strictement sous l'hypothèse d'homoscédasticité des résidus marginaux pour proposer une stratégie d'estimation qui permet d'avoir des estimations empiriquement consistantes. Cette approche que nous proposons est en fait une généralisation (extension) de la méthode d'estimation proposée par Bates et ses co-auteurs [14], et que nous avons détaillée à la Section 1.2.2. Cette méthode qui pourrait être qualifiée d'estimation à la loupe, surtout en ce qui concerne les paramètres de variance des effets aléatoires, n'est finalement qu'une réécriture habile de la vraisemblance vue comme fonction, uniquement, des paramètres de variance sur lesquels le processus d'optimisation est itéré jusqu'à convergence.

Comme expliqué un peu plus haut, nous allons nous restreindre à la dimension *k* = 2, puis-
qu'une écriture générale paraît trop lourde mais relativement simple à obtenir quand on part du
cas $k = 2$. Posons $\beta = (\beta_1^\top, \beta_2^\top)^\top$, $\sigma = (\sigma_1^2, \sigma_2^2)^\top$, $\theta = (\theta_1^\top, \theta_2^\top)^\top$, où $\sigma_1^2$ et $\sigma_2^2$ sont les variances résiduelles dans les deux dimensions respectivement, et $\theta_1$ et $\theta_2$ sont les vecteurs respectifs des composantes de variances des effets aléatoires dans les deux dimensions 1 et 2 du modèle. Notons en suite $Y_\sigma = (\sqrt{\sigma_2^2}y_1^\top, \sqrt{\sigma_1^2}y_2^\top)^\top$, $X_\sigma = \begin{pmatrix} \sqrt{\sigma_2^2}X_1 & \mathbf{0} \\ \mathbf{0} & \sqrt{\sigma_1^2}X_2 \end{pmatrix}$, $Z_{\sigma\theta} = \begin{pmatrix} \sqrt{\sigma_2^2}Z_1\Lambda_{\theta_1} & \mathbf{0} \\ \mathbf{0} & \sqrt{\sigma_1^2}Z_2\Lambda_{\theta_2} \end{pmatrix}$. Soit $u_1$ et $u_2$ les variables sphériques des effets aléatoires $\gamma_1 \in \mathbb{R}^{q_1}$ et $\gamma_2 \in \mathbb{R}^{q_2}$ de chacune des dimensions, et $q = q_1 + q_2$.

$$ \gamma_1 = \Lambda_{\theta_1} u_1, \quad \gamma_2 = \Lambda_{\theta_2} u_2 \qquad (1.74) $$

tel que

$$ u = \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} \sim \mathcal{N}(\mathbf{0}, \Sigma_u), \quad \text{avec} \quad \Sigma_u = \begin{pmatrix} \sigma_1^2 \mathbf{I}_{q_1} & \sigma_1 \sigma_2 \rho \\ \sigma_1 \sigma_2 \rho^\top & \sigma_2^2 \mathbf{I}_{q_2} \end{pmatrix}. \qquad (1.75) $$

$\rho = \text{diag}(\rho, \dots, \rho)$, où $\rho$ est la matrice qui contient les corrélations entre les effets aléatoires marginaux $\gamma_1$ et $\gamma_2$. Notons $\mu_{\mathcal{U}|\mathcal{Y}=y}$ l'espérance de $\mathcal{U}$ conditionnellement à $\mathcal{Y} = y$. La log-vraisemblance de $\beta, \sigma, \theta$ et $\rho$ conditionnellement à $y$ s'écrit

$$ \ell(\beta, \theta, \rho, \sigma | \mathbf{y}) = - \frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathcal{Y}=y}) + \|R_X(\beta - \hat{\beta}_{\theta, \rho, \sigma})\|^2}{2\sigma_1^2\sigma_2^2} - \frac{N-q}{2} \log(\sigma_1^2\sigma_2^2) \\ - \frac{1}{2} \log(|\Sigma_u|) - \frac{1}{2} \log(|L_{\theta, \rho, \sigma}|^2) \qquad (1.76) $$

où, $\hat{\beta}_{\theta, \rho, \sigma}$ et $\mu_{\mathcal{U}|\mathcal{Y}=y}$ satisfont

$$ \begin{pmatrix} X_{\sigma}^{\top} X_{\sigma} & X_{\sigma}^{\top} Z_{\sigma\theta} \\ Z_{\sigma\theta}^{\top} X_{\sigma} & Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2} \Sigma_u^{-1} \end{pmatrix} \begin{pmatrix} \hat{\beta}_{\theta, \rho, \sigma} \\ \mu_{\mathcal{U}|\mathcal{Y}=y} \end{pmatrix} = \begin{pmatrix} X_{\sigma}^{\top} \\ Z_{\sigma\theta}^{\top} \end{pmatrix} Y_{\sigma}, \qquad (1.77) $$

$$ r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathcal{Y}=y}) = \|Y_{\sigma} - X_{\sigma} \hat{\beta}_{\theta, \rho, \sigma} - Z_{\sigma\theta} \mu_{\mathcal{U}|\mathcal{Y}=y}\|^2 + \sigma_1^2 \sigma_2^2 \mu_{\mathcal{U}|\mathcal{Y}=y}^{\top} \Sigma_u^{-1} \mu_{\mathcal{U}|\mathcal{Y}=y}, \qquad (1.78) $$

$L_{\theta, \rho, \sigma}$ vérifie

$$ L_{\theta, \rho, \sigma} L_{\theta, \rho, \sigma}^{\top} = Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2} \Sigma_u^{-1}, \qquad (1.79) $$

et $R_X$ vérifie

$$ \begin{pmatrix} X_{\sigma}^{\top} X_{\sigma} & X_{\sigma}^{\top} Z_{\sigma\theta} \\ Z_{\sigma\theta}^{\top} X_{\sigma} & L_{\theta, \rho, \sigma} L_{\theta, \rho, \sigma}^{\top} \end{pmatrix} = \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta, \rho, \sigma}^{\top} \end{pmatrix}^{\top} \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta, \rho, \sigma}^{\top} \end{pmatrix}. \qquad (1.80) $$

En résolvant l'équation $\frac{\partial \ell(\beta, \theta, \rho, \sigma | y)}{\partial \beta} = 0$, on obtient $\beta = \hat{\beta}_{\theta, \rho, \sigma}$ qu'on re-injecte dans l'expression de $\ell(\beta, \theta, \rho, \sigma | y)$ pour obtenir
<table><caption>TABLE 1.3 : Erreur quadratique moyenne des estimateurs avec son intervalle de confiance à 95% estimés sur 100 replications de données simulées en prenant $n \in \{50, 60, 100, 300\}$ et $N \in \{600, 1000, 3000\}$.</caption><thead><tr><th>Parameter</th><th>n</th><th>N = 600</th><th>N = 1000</th><th>N = 3000</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50</td><td>2.43 (0.11 - 7.11)</td><td>1.89 (0.22 - 4.89)</td><td>1.02 (0.07 - 2.44)</td></tr><tr><td>60</td><td>2.57 (0.14 - 7.87)</td><td>2.13 (0.26 - 5.55)</td><td>0.77 (0.10 - 2.09)</td></tr><tr><td>100</td><td>2.27 (0.16 - 5.61)</td><td>1.55 (0.14 - 5.17)</td><td>0.71 (0.04 - 1.85)</td></tr><tr><td>300</td><td>3.16 (0.14 - 10.54)</td><td>1.70 (0.10 - 4.55)</td><td>0.51 (0.04 - 1.36)</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>50</td><td>5.50(0.09 - 16.35)</td><td>3.26 (0.02 - 11.64)</td><td>2.06 (0.09 - 5.98)</td></tr><tr><td>60</td><td>5.06 (0.12 - 15.09)</td><td>3.22 (0.02 - 10.24)</td><td>1.78 (0.10 - 5.62)</td></tr><tr><td>100</td><td>4.33 (0.03 - 13.17)</td><td>2.37 (0.02 - 6.89)</td><td>1.06 (0.02 - 3.60)</td></tr><tr><td>300</td><td>4.58 (0.18 - 15.92)</td><td>2.43 (0.05 - 7.39)</td><td>0.90 (0.04 - 2.88)</td></tr><tr><td rowspan="4">σ<sub>1</sub></td><td>50</td><td>0.03 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.09)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.11)</td><td>0.02 (0.00 - 0.08)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>100</td><td>0.03 (0.00 - 0.13)</td><td>0.01 (0.00 - 0.06)</td><td>0.00 (0.00 - 0.01)</td></tr><tr><td>300</td><td>0.05 (0.00 - 0.23)</td><td>0.03 (0.00 - 0.13)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td rowspan="4">σ<sub>2</sub></td><td>50</td><td>0.04 (0.00 - 0.15)</td><td>0.03 (0.00 - 0.08)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>100</td><td>0.06 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.11)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>300</td><td>0.13 (0.00 - 0.45)</td><td>0.04 (0.00 - 0.15)</td><td>0.01 (0.00 - 0.05)</td></tr><tr><td rowspan="4">Γ̄</td><td>50</td><td>0.90 (0.12 - 2.41)</td><td>0.62 (0.06 - 1.41)</td><td>0.45 (0.04 - 1.14)</td></tr><tr><td>60</td><td>1.07 (0.06 - 2.64)</td><td>0.68 (0.08 - 1.98)</td><td>0.25 (0.03 - 0.66)</td></tr><tr><td>100</td><td>0.90 (0.05 - 2.40)</td><td>0.47 (0.03 - 1.29)</td><td>0.21 (0.02 - 0.71)</td></tr><tr><td>300</td><td>1.45 (0.19 - 4.39)</td><td>0.69 (0.06 - 1.94)</td><td>0.17 (0.02 - 0.57)</td></tr></tbody></table>

$$ \tilde{\ell}(\theta, \rho, \sigma | y) = -\frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{u}|y=y)}{2\sigma_1^2\sigma_2^2} - \frac{N-q}{2}\log(\sigma_1^2\sigma_2^2) \\ -\frac{1}{2}\log(|\Sigma_u|) - \frac{1}{2}\log(|L_{\theta,\rho,\sigma}|^2). \qquad (1.81) $$

Les estimateurs du maximum de vraisemblance $\hat{\beta}$, $\hat{\sigma}$, $\hat{\theta}$, $\hat{\rho}$ de $\beta$, $\sigma$, $\theta$ et $\rho$ satisfont à

$$ (\hat{\theta}, \hat{\rho}, \hat{\sigma}) = \arg\max_{\theta,\rho,\sigma} \tilde{\ell}(\theta, \rho, \sigma | y) \quad \text{et} \quad \hat{\beta} = \hat{\beta}_{\hat{\theta},\hat{\rho},\hat{\sigma}} \qquad (1.82) $$

Soit $p, p_1$ et $p_2$ tels que $\dim(X_1) = N \times p_1$, $\dim(X_2) = N \times p_2$ et $p = p_1 + p_2$. Avec les mêmes notations, le critère REML à optimiser s'écrit

$$ \mathcal{L}(\sigma, \theta, \rho | y) = \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{u}|y=y)}{2\sigma_1^2\sigma_2^2}\right] (\sigma_1^2\sigma_2^2)^{\frac{p+q-N}{2}}}{(2\pi)^{(2N-p)/2} |\Sigma_u|^{1/2} |L_{\theta,\rho,\sigma}| |R_X|} \qquad (1.83) $$

Sur 100 jeux de données simulées, nous avons investigué empiriquement les performances de cette stratégie d'estimation en calculant l'erreur quadratique moyenne de chaque estimateur sur les 100 réplications. Les résultats obtenus, contenus dans le tableau 1.3 montrent clairement une nette amélioration de cette erreur à tous les niveaux, comparativement aux EM-estimateurs. Quand on regarde particulièrement les paramètres de variance des effets aléatoires, on constate que la présente méthode est de loin plus précise que l'algorithme EM. Avec des jeux de données
<table><caption>TABLE 1.4 : Erreur quadratique moyenne des estimateurs avec son intervalle de confiance à 95% estimés sur 100 replications de données simulées en prenant (n, N) ∈ {(500, 7000), (600, 8000), (800, 10000), (1000, 15000)}.</caption><thead><tr><th>Parameter</th><th>n = 500, N = 7000</th><th>n = 600, N = 8000</th><th>n = 800, N = 10000</th><th>n = 1000, N = 15000</th></tr></thead><tbody><tr><td>β<sub>1</sub></td><td>0.22 (0.01 - 0.62)</td><td>0.17 (0.01 - 0.48)</td><td>0.16 (0.01 - 0.47)</td><td>0.11 (0.00 - 0.32)</td></tr><tr><td>β<sub>2</sub></td><td>0.32 (0.00 - 1.01)</td><td>0.34 (0.01 - 1.16)</td><td>0.25 (0.02 - 0.67)</td><td>0.21 (0.00 - 0.69)</td></tr><tr><td>σ<sub>1</sub></td><td>0.00 (0.00 - 0.00)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.00)</td><td>0.00 (0.00 - 0.00)</td></tr><tr><td>σ<sub>2</sub></td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.00)</td></tr><tr><td>Γ</td><td>0.09 (0.00 - 0.25)</td><td>0.07 (0.00 - 0.19)</td><td>0.06 (0.00 - 0.19)</td><td>0.03 (0.00 - 0.09)</td></tr></tbody></table>

de plus grandes tailles (voir tableau 1.4), il devient très clair que cette méthode permet d'obtenir des estimateurs asymptotiquement consistants.

## Contributions relatives au chapitre 4 : Procédure de sélection d'effets fixes dans le modèle linéaire à effets mixtes : Approcher les performances d'une pénalité de type $L_0$

La procédure de sélection des effets fixes que nous proposons dans ce chapitre se base sur l'écriture de la vraisemblance proposée par Bates et ses co-auteurs [14] avec une légère modification, en présence d'un terme de pénalité de type *adaptive ridge*. L'exposé que nous faisons ici concerne seulement les modifications apportées à la méthode, quant à la façon d'écrire la vraisemblance. Puisqu'il est question de faire de la sélection parmi les effets fixes, nous avons réécrit la vraisemblance comme fonction à la fois des effets fixes et des paramètres de variance ; ce qui n'est pas le cas de l'approche initiale où la vraisemblance apparaît comme fonction uniquement des paramètres de variance des effets aléatoires, en présence d'une mise à jour itérative de "l'expression" des effets fixes et des effets aléatoires. Avec les mêmes notations qu'à la Section 1.2.2, on a

$$ \|y - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2 = \left\| \begin{pmatrix} y - X\beta \\ \mathbf{0} \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} u \right\|^2 \quad (1.84) $$

$$ \mu_{u|y=y} = \arg\min_{u} \left\| \begin{pmatrix} y - X\beta \\ \mathbf{0} \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} u \right\|^2 \\ \implies \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix}^{\top} \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} \mu_{u|y=y} = \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix}^{\top} \begin{pmatrix} y - X\beta \\ \mathbf{0} \end{pmatrix}, \quad (1.85) $$

et $\|y - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2$ se réécrit comme suit

$$ \|y - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2 = \|y - X\beta - Z\Lambda_{\theta}\mu_{u|y=y}\|^2 + \|\mu_{u|y=y}\|^2 + $$
$$
\begin{align}
& (u - \mu_{U|Y=y})^{\top} \underbrace{( (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + I_q )}_{=L_{\theta}^{\top} L_{\theta}} (u - \mu_{U|Y=y}) \nonumber \\
& = \|y - X\beta - Z\Lambda_{\theta}\mu_{U|Y=y}\|^2 + \|\mu_{U|Y=y}\|^2 + \|L_{\theta}(u - \mu_{U|Y=y})\|^2. \tag{1.86}
\end{align}
$$

Ce qui fait que la vraisemblance s'écrit

$$
L(\beta, \sigma^2, \theta|y) = |L_{\theta}|^{-1} (2\pi\sigma^2)^{-n/2} \exp \left[ - \frac{\|y - X\beta - Z\Lambda_{\theta}\mu_{U|Y=y}\|^2 + \|\mu_{U|Y=y}\|^2}{2\sigma^2} \right]. \quad (1.87)
$$

Le critère à optimiser pour la sélection d'effets fixes s'écrit

$$
C(\beta, \sigma^2, \theta|y) = \log (L(\beta, \sigma^2, \theta|y)) - \frac{\lambda}{2} \beta^\top W \beta, \quad (1.88)
$$

où, W est une matrice de poids à mettre à jour au cours de la procédure de sélection. $\frac{\partial C(\beta, \sigma^2, \theta|y)}{\partial \sigma^2} = 0 \implies \hat{\sigma}^2 = \frac{\|y - X\beta - Z\Lambda_\theta \mu_{U|Y=y}\|^2 + \|\mu_{U|Y=y}\|^2}{n}$. Cette expression de $\hat{\sigma}^2$ re-injectée dans celle de $C(\beta, \sigma^2, \theta|y)$ nous fournit

$$
-2\tilde{C}(\beta, \theta|y) = \log|L_{\theta}|^2 + n \left[ 1 + \log \left( \frac{2\pi(\|y - X\beta - Z\Lambda_{\theta}\mu_{U|Y=y}\|^2 + \|\mu_{U|Y=y}\|^2)}{n} \right) \right] + \lambda\beta^{\top}W\beta \quad (1.89)
$$

qui est le critère qui sera minimisé pour faire la sélection.
# MULTIVARIATE LONGITUDINAL ANALYSIS WITH BIVARIATE CORRELATION TEST

Dans le cadre de l'ajustement d'un modèle linéaire à effets mixtes à des données longitudinales multidimensionnelles, nous proposons ici des expressions plus générales (que celles rencontrées dans la littérature) des estimateurs du maximum de vraisemblance des paramètres, utilisant l'algorithme EM. Ces estimateurs sont aussi directement utiles dans le cadre plus général de l'analyse de données multidimensionnelles multi-niveaux. L'idée générale de l'analyse jointe de plusieurs variables statistiques par un modèle multidimensionnel à effets mixtes, permettant de modéliser la distribution jointe de tous les effets aléatoires spécifiques à chaque dimension du modèle, est toujours applicable quel que soit le nombre de variables dépendantes. Cependant, la dimension du vecteur de tous les effets aléatoires du modèle peut très vite croître avec le nombre de variables à l'analyse, entrainant ainsi de sérieux problèmes numériques d'optimisation. Dans ce chapitre qui a fait l'objet d'une publication dans PloS One, nous proposons de tester la dépendance entre deux dimensions du modèle par un test de rapport de vraisemblance qui porte sur la significativité globale des corrélations entre les effets aléatoires de ces dimensions. Cette méthode nous permettra de construire un modèle plus parcimonieux en terme de nombre de paramètres de variance relatifs aux effets aléatoires, à travers une procédure de sélection pas-à-pas ascendante. Ce test de corrélation bi-dimensionnel construit ici a montré ses performances dans les études de simulations que nous avons menées, mais aussi son utilité dans l'analyse de données réelles.
# Multivariate Longitudinal Analysis with Bivariate Correlation Test

Eric Houngla Adjakossa¹,²*, Ibrahim Sadissou³,⁴, Mahouton Norbert Hounkonnou², Gregory Nuel¹

1 Laboratoire de Probabilités et Modèles Aléatoires /Université Pierre et Marie Curie, Case courrier 188 - 4, Place Jussieu 75252 Paris cedex 05 France

2 University of Abomey-Calavi, 072 B.P. 50 Cotonou, Republic of Benin

3 Laboratoire de Biologie et de Physiologie Cellulaires /University of Abomey-Calavi, Cotonou, Republic of Benin

4 Centre d'Etude et de Recherche sur le Paludisme Associé à la Grossesse et à l'Enfance (CERPAGE), Cotonou, Republic of Benin

* ericadjakossah@gmail.com

## Abstract

In the context of multivariate multilevel data analysis, this paper focuses on the multivariate linear mixed-effects model, including all the correlations between the random effects when the dimensional residual terms are assumed uncorrelated. Using the EM algorithm, we suggest more general expressions of the model's parameters estimators. These estimators can be used in the framework of the multivariate longitudinal data analysis as well as in the more general context of the analysis of multivariate multilevel data. By using a likelihood ratio test, we test the significance of the correlations between the random effects of two dependent variables of the model, in order to investigate whether or not it is useful to model these dependent variables jointly. Simulation studies are done to assess both the parameter recovery performance of the EM estimators and the power of the test. Using two empirical data sets which are of longitudinal multivariate type and multivariate multilevel type, respectively, the usefulness of the test is illustrated.

## 2.1 Introduction

In statistical studies, one often needs to analyze data with nested sources of variability: e.g., pupils in classes, employees in companies, repeated measurements in subjects, etc. [140] referred to these type of data as grouped data which are also named multilevel data, hierarchical data or nested data in the literature [72, 177, 210]. In the analysis of such data, it is usually illuminating to take account of the variability associated with each level of nesting. There is variability, e.g., between pupils but also between classes. The measurements related to a specific subject (level of nesting) can be correlated, while observations from different subjects are usually independent, and one may draw wrong conclusions if either of these sources of variability is ignored [206]. A series of works
in statistical literature focus on the analysis of univariate multilevel data (or univariate grouped data) where a single outcome of interest is analyzed [46, 110, 114, 129, 193, 204]. Such analyses are generally simple to deal with due to the availability of many software packages conceived to perform them [16, 117, 142]. In practice, many scientific questions of interest require to focus on multiple outcomes, all arising from the same multilevel study, leading to the so-called multi-variate multilevel data. For example, to answer some questions of interest, [57] analyzed hearing threshold data (in the Baltimore Longitudinal Study on Aging) [174] which consisted in the longitudinal recording of 22 variables. [185] also studied the joint evolution of HIV RNA and CD4+ T lymphocytes in a cohort of HIV-1 infected patients treated with highly active antiretroviral treatment, by jointly analyzing both markers. [182] used multivariate multilevel regression analysis to investigate individual level determinants of self rated health and happiness, as well as the extent of community level covariation in health and happiness. [190] also used multivariate multilevel analysis to jointly model three commonly used indicators of fear of crime which are: feeling unsafe alone at home after dark, feeling unsafe walking alone after dark and worry about becoming a victim of crime. A variety of works were devoted to joint modeling during the last few decades (see e.g., [26, 32, 59, 183, 198]).

These analyses often require a specification of the joint density of all outcomes or, at least, the correlation structure of the data and therefore can lead to the parsimony and/or computation (optimization) problems as well as to numerical difficulties in statistical inference, when the dimension of these outcomes increases. Many analysis strategies were proposed in the statistical literature to circumvent these problems. These strategies generally consist in reducing the dimensionality of the multivariate vector of outcomes and/or in using a small number of latent variables to model correlations within these data. Joint analysis of multivariate multilevel data then requires a trade-off between the increase of the computational complexity and the gain in information.

In this work, we focus on the multivariate linear mixed-effects model, including all the correlations between the random effects along with the independent marginal (dimensional) residuals. The correlations between two dependent variables are then those from the random effects related to these dependent variables. The class of mixed-effects models considered here assumes that both the random effects and the errors (residuals) follow Gaussian distributions. These models are intended for the analysis of multivariate multilevel data in which the dependent variables are continuous.

We use the EM algorithm to estimate the parameters of the model but here, we have two novelties: 1) we suggest a general expression of EM-based estimators which can help in analyzing multivariate longitudinal data as well as the multivariate multilevel data, not of the longitudinal type, and 2) we test the significance of the correlations between the random effects of two dependent variables, using the likelihood ratio test which allows to decide if some dependent variables are significantly correlated or not. By using this bivariate correlation test, the novelty here is the illustration, through empirical data, of some of the consequences of performing separate analyses when a joint analysis is required. Two dependent variables which are found to be uncorrelated after this test will be analyzed with two independent models (or analyzed separately). This strategy may be considered as a way toward the obtaining of a more parsimonious model in high dimension without losing much information. It may also be used in a joint modeling selection procedure.
The paper is organized as follows. In Section 2, contributions of previous works are briefly presented. We also present in this section the EM-based estimators of the parameters of the multivariate linear mixed model. Simulations studies are done in Section 3 where we also discuss the power of the likelihood ratio test which allows to test the significance of the correlation between two response variables. Two illustrations on empirical data are also done in Section 3. The first, concerning bivariate two-level data, is about a study on the effects of school differences on pupils' progress in Dutch language and arithmetics in the Netherlands. The second illustration concerns a longitudinal study on the immune response to malaria of infants in Benin.

## 2.2 Materials and Methods

### 2.2.1 Previous works

In this Section, we briefly recall the framework of the multivariate multilevel analysis (see for instance, [10, 194]). We can basically distinguish two main approaches to model such data: those which specify the joint distribution of all outcomes without the use of latent structures, and the models using latent structures. We denote by $y_1, \dots, y_m$ the $m$ dependent vectors of interest, and $\mathbf{y} = (y_1^\top, \dots, y_m^\top)^\top$.

**Modeling methods without latent structures**

The first approach which is that of the modeling without latent structures comprises three sub- approaches consisting in a) direct specification of the correlation structure of **y**, b) analysis without explicit modeling of the correlation structure of **y** and c) conditional models.

In the case of direct specification of Cov(y), [68] and [135] factorized the covariance matrix of y by using the Kronecker product in order to have more parsimonious models in the context of fully balanced data. With the same idea of having a parsimonious structure, [30] specified the intra-outcome and inter-outcome correlations, respectively, as follows: Corr(y_k(t), y_k(s)) = exp(α|s - t|θ_k) and Corr(y_k(t), y_k'(s)) = exp((α|s - t| + 1)θ_kk'), with t and s indicating the time, and k and k' indicating the dimension. Although these models are useful, they are often too restrictive and may not be realistic in many applications, especially when the data, for example, in the longitudinal studies are unbalanced (i.e. the number of available measurements per subject and the time points at which the measurements were taken often differ from one subject to another). Another class of joint models, specifying directly the joint distribution of y, and whose application is often not straightforward, due to unbalanced data structures is the so-called copula model [132, 175]. Denoting by $F_i, i = 1, \dots, m$ the cumulative distribution function of the ith component of y, $y_i$, a copula model is defined by an $m$-dimensional cumulative distribution function $C(u_1, \dots, u_m)$ with uniform marginals such that $F(y_1, \dots, y_m) = C(u_1, \dots, u_m)$ with $(U_1, \dots, U_m) = (F_1(y_1), \dots, F_m(y_m))$, where F is the joint cumulative distribution of y. That is, the joint distribution of y can be written in terms of its marginal distributions and a copula which describes the dependence structure between its components. While the construction of copulas is mathematically elegant, parameters estimation is often not feasible, especially in high-dimensional situations [194]. One of the rare applications of the copula-based modeling in the multivariate mul-
tilevel data analysis framework was proposed by [106], who studied the hemodynamic effect of a new antidepressant on the diastolic blood pressure, the systolic blood pressure and the heart rate of 10 healthy volunteers. They separately modeled, at first, each longitudinal series of response and used a copula to relate the marginal distributions of these responses at each observation time. In a second step, at each observation time, the conditional (on the past) distributions of each response were related using another copula describing the relationship between the corresponding variables. One of the advantages of this approach is that there is no need to use the same family of distributions for all response variables. As [121] used ARIMA process to model the error structure of earnings in a longitudinal data analysis context, time series models can also be used for modeling multivariate multilevel data in order to describe the dynamic dependence between variables and perform forecasting. The most commonly used multivariate time series model, the vector autoregressive (VAR) model which is relatively easy to estimate, is found to be similar to the multivariate multiple linear regression [187] where the errors for different response variables on the same trial are set to be correlated [99]. Other examples of VAR modeling include [188] and [189], but one drawback of the model is that the number of parameters can become very large, potentially leading to estimation problems [93].

Regarding analysis without explicit specification of Cov(y), [110] proposed an extension of generalized linear models to the analysis of longitudinal data, where they introduced a class of estimating equations called generalized estimating equations (GEE). GEE estimation ensures consistent estimates of the regression parameters without specifying the joint distribution of a subject's observations. That is, GEE replaces $V[y]$ by the so-called working covariance matrix $W(\alpha)$ which depends on an unknown vector $\alpha$ to estimate. The related working correlation matrix, $R(\alpha)$, is also considered. Incorrect choice of $W(\alpha)$ does not affect the consistency of the regression parameters' estimators [110]. [111] discussed the use of GEEs with multivariate discrete variables, where focus was on the modeling of the marginal (dimensional) means of these variables and their pairwise associations. The extension of the GEE method to mixed continuous-discrete responses was discussed by [205] and [143]. [156] also avoided the need of explicit modeling of the covariance structure of bivariate longitudinal responses by using SUR [206] and GEE. As pointed out by [40], ambiguities concerning the definition of the working covariance matrix can result in a breakdown of the GEE-based estimation. For example in the longitudinal data analysis, if the true structure of correlation is equicorrelation, $(R_i^*)_{jk} = \rho$, and that the working structure is autoregressive, $(R_i)_{jk} = \alpha^{|j-k|}$, there is no solution for $\hat{\alpha}$ when $-1/2 \le \rho < -1/3$ [40]. This can be viewed as the major drawback of the GEE method since it can lead to the misspecification of within-subject associations in the context of longitudinal data analysis, for instance. Examples of procedures which bypass the need to explicitly model the underlying covariance structure of y include [79, 80, 156]. These procedures, generally, consist in regressing each component of y on relevant covariates of interest, followed by combination of these regression coefficients into a single global estimate of the covariates effect [10].

One way to avoid the direct specification of the joint distribution of y is to factorize it, leading to the so-called conditional models [73]. For two responses, the joint density $f(y_1, y_2)$ can be written as follows:

$$f(y_1, y_2) = f(y_1|y_2)f(y_2) = f(y_2|y_1)f(y_1) \qquad (2.1)$$
The choice of the conditioning response is of course arbitrary and requires very careful reflection about plausible associations between components of $y$. For example, in the specification of a conditional model such as $f(y_1|y_2)$, $y_2$ plays the role of covariate and different choices can lead to completely opposite results and conclusions [208]. In a clinical trial, for example, none of these factorizations will be of interest due to the conditioning on a post-randomization outcome which may partially attenuate the treatment effect on the other [194]. Another drawback of conditional models is that they do not directly lead to marginal inferences. Suppose that scientific interest would be in a comparison of the rate of longitudinal change in average of $y_1$ and $y_2$. The factorization $f(y_1, y_2) = f(y_1|y_2)f(y_2)$ directly allows for inferences about the marginal evolution of $y_2$, but the marginal expectation of $y_1$ requires computation of $\mathbb{E}[y_1] = \mathbb{E}[\mathbb{E}[y_1|y_2]]$, which, depending on the actual models, may be far from straightforward [194].

## Modeling methods using latent structures

The second approach regarding models using latent structures can also be split in two sub-approaches including the strategy based on the reduction of the dimensionality of $y$ and the mixed-effect models. The general idea of reducing the dimension of $y$ is to use principal-component type analysis, or a summary function, to first reduce the dimensionality of $y$ and then, use standard univariate multilevel models for the analysis of the principal factors or the retained summaries of $y$ [49, 58, 84, 123, 137]. Although it is useful, simple to understand and easy to compute, this strategy of dimension reduction has some drawbacks such as the loss of information as discussed by [10] and [194]. [10] used this approach and retained the first principal-component only which explains 31% of the total variation in their data. They found out that the summary function does not have any physical significance and the inference results cannot be interpreted in terms of the effect of the covariates on the original (response) variables. They also found that the method fails to explore the association of the components of $y$ along time, in the case of longitudinal studies. Furthermore, the method is not applicable in situations where all the components of $y$ are not measured at the same time point [10], although a possible extension might be the use of functional principal components [148].

Regarding the mixed-effect models, [57], [151], [120], [152] and [18] proposed the use of random-effects models for multivariate longitudinal data. They pointed out that the main disadvantage of joining separate mixed models by allowing their model-specific random effects to be correlated is the increase of the dimension of the total vector of random effects with the number of outcomes, leading to computational problems. To circumvent these problems, [57] noted that all parameters in the joint model can be estimated by fitting all the bivariate models, based on

$$f(y_s, y_t) = \iint f(y_s|\gamma_s)f(y_t|\gamma_t)f(\gamma_s, \gamma_t)d\gamma_s d\gamma_t$$

for all $m(m-1)/2$ pairs $(y_s, y_t)$, $1 \le s < t \le m$, resulting from the main multivariate model. Estimators for the main parameters are obtained by averaging over the results from fitting the $m(m-1)/2$ pairwise models. They then showed that the pseudo-likelihood theory can be used to derive the asymptotic distribution of these estimators, and used SAS procedures for mixed models [117] based on the Newton-Raphson algorithm to fit their models, following the approach in [185].
In some multilevel studies, focus is not to directly model $\mathbf{y}$, but a few number of latent variables which cannot be quantified directly (e.g., depression and anxiety), but through measurements of $\mathbf{y}$. In such situations, analysis may be conducted in two steps: the first produces the obtaining of the latent variables and the second proceeds to the joint analysis of these latent variables. For example, [6] proposed a latent factor linear mixed model to capture the joint trend over time of latent variables. the authors reduced, indeed, the high-dimensional responses to low-dimensional latent factors by the factor analysis model, and then used the multivariate linear mixed model to study the longitudinal trends of these latent factors, where the estimates have been done using the EM algorithm. To deal with missing values in multivariate longitudinal analysis using multivariate linear mixed-effects model, [161] proposed multiple imputations using Markov chain Monte Carlo, where they used EM algorithm for the parameters estimation. Here, the authors sped up the EM algorithm by analytically integrating the random effects out of the likelihood function, avoiding to treat them as missing data. [172] used EM based modeling to estimate the parameters of the multivariate linear mixed model under a SAS macro program encoded in IML.

Although the EM algorithm is known to be slow, one of the biggest advantages of this method is that it is not computationally expensive, even with a large number of response variables. In this context, our contribution is the writing of the EM-based estimators in a more general form than those used in [6], [161] and [172]. The expressions of the EM-based estimators used in this paper can easily perform any analysis in the framework of the multivariate multilevel data analysis using multivariate linear mixed-effects model.

Another technique somewhat close to those discussed in [6] is the structural equations-based techniques. For example, [19] developed linear structural equations with latent variables approach. Considering $\mathbf{y} = (y_1^\top, y_2^\top)^\top$, this approach can be expressed as follows: $y_i = \mu_i + G_i\eta_i; i = 1, 2$ and $\beta\eta_1 = \gamma\eta_2$, where $\eta_i, i = 1, 2$ are the latent variables, $\beta(m \times m)$ and $\gamma(m \times n)$ are coefficient matrices governing the linear relations of all variables involved in the $m$ structural equations. $G_i, i = 1, 2$ are known matrices. The parameters of the model may be estimated by gradient and quasi-Newton methods, or a Gauss-Newton algorithm that obtains least-squares, generalized least-squares, or maximum likelihood estimates. One modeling strategy which fuses together mixed-effects model and VAR model in order to analyze multivariate multilevel data is the so-called multilevel-VAR method. For example, [25] used the multilevel-VAR model in the context of network inference in psychopathology, where they used the population standard deviation of the person-specific random effects to construct a network representing individual variability. Examples of multilevel-VAR modeling include [67] and [93].

State space models [83] which are useful to investigate the dynamical properties of latent variables can also be used to analyze multivariate multilevel data. For example, [119] introduced an extension of the basic state space model which is flexible and general in the sense of it is applicable to any time series for multiple systems.

Methods for estimating the connectivity maps containing heterogeneity may also be applied to analyze multivariate multilevel data. [70] presented the Group Iterative Multiple Model Estimation (GIMME) approach, which addresses the issue of heterogeneity (the need for individual-level maps) in effective connectivity mapping while capitalizing on shared information to arrive at group inferences. Unlike mixed-effects models, GIMME allows for the structure of the connectivity maps to be unique across individuals [70].
One can also use a nonparametric function $f$ to handle the relationship between the components of $y$ and the covariates [55, 153, 201]. This strategy requires also to have sufficient data per subject, in the case of multivariate longitudinal data. Other estimation strategies implemented under softwares and discussed by [42] can perhaps be extended to the multivariate analysis case, when necessary.

Let us finally point out that the software packages which can easily and accurately analyze (jointly) the data of multivariate multilevel type are extremely rare, and one arranges the data and manipulates packages primarily designed for fitting univariate models to handle their analysis. The SabreR [39] package, under the R software [145], which has been devoted to jointly fitting up to three mixed-effects models, with random intercepts only, has been recently removed from the depot. These facts prove by themselves that the analysis of multivariate multilevel data in a single framework is a challenging task. Bayesian-based approaches can be implemented using packages like R2WinBUGS [181] under the R software, and are useful but very time consuming and require a good expertise from the user who can easily be discouraged.

### 2.2.2 Model and notations

The model discussed here is the multivariate linear mixed-effects model (or the multivariate linear multilevel model), including all the correlations between the random effects, but the marginal residual terms are assumed to be uncorrelated. For a more general multivariate linear mixed-effects model, the dependent variables are assumed to be correlated, conditional on the random-effects. That is, the marginal residual terms are correlated. In this paper, as in many other works (see for example, [172], [161], [56] and [6]), we assume that conditional on the random-effects, the dependent variables are uncorrelated. In the context of using EM algorithm in estimating the model parameters, this assumption allows to derive the EM-based estimators for the residual variance parameters. If the dimensional residual terms are assumed to be correlated, the EM-based estimators of theirs variance parameters are not easy to deal with and we don't treat this case here. This model assumes that both the random effects and the residuals follow Gaussian distribution, and is intended for the analysis of multivariate multilevel data in which the dependent variables are continuous. For the sake of simplicity we focus on the bivariate case ($m = 2$) in most of the paper, but the generalization to higher dimensions ($m > 2$) is straightforward. The model is as follows:

$$ 
\begin{aligned}
y_1 &= X_1\beta_1 + Z_1\gamma_1 + \varepsilon_1 \\
y_2 &= X_2\beta_2 + Z_2\gamma_2 + \varepsilon_2
\end{aligned} 
\qquad (2.2) 
$$

$$ \gamma = \begin{pmatrix} \gamma_1 \\ \gamma_2 \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \mathbf{\Gamma} = \begin{pmatrix} \Gamma_1 & \Gamma_{12} \\ \Gamma_{12}^\top & \Gamma_2 \end{pmatrix} \right), \qquad (2.3) $$

$$ \varepsilon = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \mathbf{\Sigma} = \begin{pmatrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{pmatrix} \right); \quad \gamma \perp \varepsilon \qquad (2.4) $$
For $k \in \{1, 2\}$, $\beta_k$ and $\gamma_k$ denote respectively the fixed effects and the random effects vector of covariates, while $\varepsilon_k$ is the residual component. $X_k$ is a matrix of covariates and $Z_k$ a covariates-based design matrix. $\dim(X_k)=N_k \times p_k$ and $\dim(Z_k)=N_k \times q_k$, where $N_k$ is the total number of observations in the dimension $k$ of the model. $p_k$ and $q_k$ are, respectively, the number of fixed effect related covariates and the number of random effect related covariates in the dimension $k$ of the model. If $N_k$ is a constant $N$ for any $k$, the index $k$ will be removed and $N$ will denote the total number of observations in all dimensions of the model. The bold symbols represent parameters of multiple dimensions (i.e. $\Sigma_1$ concerns dimension 1 of the model while $\Sigma$ concerns both dimensions).

Another way to easily understand the model is to express it using the levels of the covariate related to the random-effects. This expression (subject-based version) of the model is, generally, used in the framework of longitudinal data analysis, and lead to EM-based estimators (expressions) which are a particular case of the estimators expressions obtained in Equations 2.17, 2.18 and 2.19 (for example, see [172]). Denoting by $n$ the total number of subjects involved in the longitudinal study, the model can be expressed as follows:

denoting by $i$ a subject, for $i = 1, \dots, n$

$$ 
\begin{aligned}
y_{1i} &= X_{1i}\beta_1 + Z_{1i}\gamma_{1i} + \varepsilon_{1i} \\
y_{2i} &= X_{2i}\beta_2 + Z_{2i}\gamma_{2i} + \varepsilon_{2i}
\end{aligned}
\qquad (2.5) $$

with

$$ \gamma_i = \begin{pmatrix} \gamma_{1i} \\ \gamma_{2i} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \bar{\mathbf{\Gamma}} = \begin{pmatrix} \bar{\Gamma}_1 & \bar{\Gamma}_{12} \\ \bar{\Gamma}_{12}^\top & \bar{\Gamma}_2 \end{pmatrix} \right) \qquad (2.6) $$

and

$$ \varepsilon_i = \begin{pmatrix} \varepsilon_{1i} \\ \varepsilon_{2i} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}; \bar{\mathbf{\Sigma}} = \begin{pmatrix} \sigma_1^2 \mathbf{I}_{N_{1i}} & 0 \\ 0 & \sigma_2^2 \mathbf{I}_{N_{2i}} \end{pmatrix} \right) \qquad (2.7) $$

$N_{1i}$ and $N_{2i}$ are the dimensions of $y_{1i}$ and $y_{2i}$, respectively. Here, we assume that the marginal residuals are homoscedastic ($\mathbb{V}[\varepsilon_{ki}] = \sigma_k^2 I_{N_{ki}}$, $k = 1, 2$), but the residual covariance matrices can be of full form as in Equation 2.4. In order to make clear the relation between the model described by Equations 2.2, 2.3 and 2.4, and its version expressed by Equations 2.5, 2.6 and 2.7, we propose below a detailed example.

**Detailed example** We place ourselves in the case of longitudinal data where we observe two response variables $y_1$ and $y_2$ which are respectively the weight (kg) and the size (cm) of infants according to the score ($V_2$) of the quality of their food as well as the quality score ($V_1$) of their mothers' food. Infants are $n = 3$ girls (sex = F) and boys (sex = M) who are monitored over time. The dataset is presented by Table 2.1.

Suppose that the model at each of two dimensions has one random intercept by subject (infant) and one random slope by subject in the direction of the infant's age (in months). For example, considering an identifiability constraint covering the sex variable whose level F is the reference, the bivariate linear mixed model can be written as follows:
<table><caption>Table 2.1: Example of data</caption><thead><tr><th rowspan="2">subject</th><th rowspan="2">age</th><th rowspan="2">sex</th><th colspan="4">Response variables</th></tr><tr><th>V<sub>1</sub></th><th>V<sub>2</sub></th><th>y<sub>1</sub></th><th>y<sub>2</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>F</td><td>63.76</td><td>38.16</td><td>3.14</td><td>47.82</td></tr><tr><td>2</td><td>4</td><td>M</td><td>100.88</td><td>41.46</td><td>4.87</td><td>64.02</td></tr><tr><td>1</td><td>6</td><td>F</td><td>60.98</td><td>41.37</td><td>8.43</td><td>73.21</td></tr><tr><td>3</td><td>0</td><td>M</td><td>93.24</td><td>48.76</td><td>2.82</td><td>44.93</td></tr><tr><td>2</td><td>7</td><td>M</td><td>101.95</td><td>44.79</td><td>8.03</td><td>89.54</td></tr><tr><td>2</td><td>10</td><td>M</td><td>99.24</td><td>48.17</td><td>10.08</td><td>92.14</td></tr><tr><td>1</td><td>16</td><td>F</td><td>NA</td><td>44.79</td><td>13.96</td><td>86.12</td></tr><tr><td>3</td><td>9</td><td>M</td><td>88.38</td><td>47.91</td><td>8.47</td><td>86.42</td></tr></tbody></table>

$$ 
\begin{aligned}
y_1 &= \underbrace{(1 \ V_1 \ sex=M)}_{=X_1} \beta_1 + Z_1 \gamma_1 + \varepsilon_1 \\
y_2 &= \underbrace{(1 \ V_2 \ sex=M)}_{=X_2} \beta_2 + Z_2 \gamma_2 + \varepsilon_2
\end{aligned} 
\qquad (2.8) 
$$

where, explicitly

$$ X_1 = \begin{pmatrix} 1 & 63.76 & 0 \\ 1 & 100.88 & 1 \\ 1 & 60.98 & 0 \\ 1 & 93.24 & 1 \\ 1 & 101.95 & 1 \\ 1 & 99.24 & 1 \\ 1 & 88.38 & 1 \end{pmatrix}; \quad X_2 = \begin{pmatrix} 1 & 38.16 & 0 \\ 1 & 41.46 & 1 \\ 1 & 41.37 & 0 \\ 1 & 48.76 & 1 \\ 1 & 44.79 & 1 \\ 1 & 48.17 & 1 \\ 1 & 44.79 & 0 \\ 1 & 47.91 & 1 \end{pmatrix} \qquad (2.9) $$

$$ Z_1 = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 4 & 0 & 0 \\ 1 & 6 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 7 & 0 & 0 \\ 0 & 0 & 1 & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 9 \end{pmatrix}; \quad Z_2 = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 4 & 0 & 0 \\ 1 & 6 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 7 & 0 & 0 \\ 0 & 0 & 1 & 10 & 0 & 0 \\ 1 & 16 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 9 \end{pmatrix} \qquad (2.10) $$

$$ \dim(X_1) = 7 \times 3 \text{ and } \dim(X_2) = 8 \times 3; \quad \dim(Z_1) = 7 \times 6 \text{ and } \dim(Z_2) = 8 \times 6 $$

In the present example we have $\dim(X_1) \neq \dim(X_2)$ and $\dim(Z_1) \neq \dim(Z_2)$ due to the presence of the NA (Not Available) within the values of the variable $V_1$. Removing information
related to this NA in the dimension 1 of the model does not affect its dimension 2.

$$ \Gamma = \begin{pmatrix} \Gamma_1 & \Gamma_{12} \\ \Gamma_{12}^\top & \Gamma_2 \end{pmatrix} $$

with,

$$ \Gamma_1 = \begin{pmatrix} \eta_1^2 & \rho_\eta \eta_1 \eta_2 & \cdot & \cdot & \cdot & \cdot \\ \rho_\eta \eta_1 \eta_2 & \eta_2^2 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \eta_1^2 & \rho_\eta \eta_1 \eta_2 & \cdot & \cdot \\ \cdot & \cdot & \rho_\eta \eta_1 \eta_2 & \eta_2^2 & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \eta_1^2 & \rho_\eta \eta_1 \eta_2 \\ \cdot & \cdot & \cdot & \cdot & \rho_\eta \eta_1 \eta_2 & \eta_2^2 \end{pmatrix} \quad (2.11) $$

$$ \Gamma_2 = \begin{pmatrix} \tau_1^2 & \rho_\tau \tau_1 \tau_2 & \cdot & \cdot & \cdot & \cdot \\ \rho_\tau \tau_1 \tau_2 & \tau_2^2 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \tau_1^2 & \rho_\tau \tau_1 \tau_2 & \cdot & \cdot \\ \cdot & \cdot & \rho_\tau \tau_1 \tau_2 & \tau_2^2 & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \tau_1^2 & \rho_\tau \tau_1 \tau_2 \\ \cdot & \cdot & \cdot & \cdot & \rho_\tau \tau_1 \tau_2 & \tau_2^2 \end{pmatrix} \quad (2.12) $$

$$ \Gamma_{12} = \begin{pmatrix} \rho_{\eta_1 \tau_1} \eta_1 \tau_1 & \rho_{\eta_1 \tau_2} \eta_1 \tau_2 & \cdot & \cdot & \cdot & \cdot \\ \rho_{\eta_2 \tau_1} \eta_2 \tau_1 & \rho_{\eta_2 \tau_2} \eta_2 \tau_2 & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \rho_{\eta_1 \tau_1} \eta_1 \tau_1 & \rho_{\eta_1 \tau_2} \eta_1 \tau_2 & \cdot & \cdot \\ \cdot & \cdot & \rho_{\eta_2 \tau_1} \eta_2 \tau_1 & \rho_{\eta_2 \tau_2} \eta_2 \tau_2 & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \rho_{\eta_1 \tau_1} \eta_1 \tau_1 & \rho_{\eta_1 \tau_2} \eta_1 \tau_2 \\ \cdot & \cdot & \cdot & \cdot & \rho_{\eta_2 \tau_1} \eta_2 \tau_1 & \rho_{\eta_2 \tau_2} \eta_2 \tau_2 \end{pmatrix} \quad (2.13) $$

$\rho_\eta, \rho_\tau, \rho_{\eta_1 \tau_1}, \rho_{\eta_1 \tau_2}, \rho_{\eta_2 \tau_1}$, and $\rho_{\eta_2 \tau_2}$ lie in $[-1, 1]$. All other parameters involved in $\Gamma_1, \Gamma_2$ and $\Gamma_{12}$ are positive real numbers.

Referring to the subject-based version of the model,

$$ \bar{\Gamma}_1 = \begin{pmatrix} \eta_1^2 & \rho_\eta \eta_1 \eta_2 \\ \rho_\eta \eta_1 \eta_2 & \eta_2^2 \end{pmatrix}, \quad \bar{\Gamma}_2 = \begin{pmatrix} \tau_1^2 & \rho_\tau \tau_1 \tau_2 \\ \rho_\tau \tau_1 \tau_2 & \tau_2^2 \end{pmatrix}, \quad \bar{\Gamma}_{12} = \begin{pmatrix} \rho_{\eta_1 \tau_1} \eta_1 \tau_1 & \rho_{\eta_1 \tau_2} \eta_1 \tau_2 \\ \rho_{\eta_2 \tau_1} \eta_2 \tau_1 & \rho_{\eta_2 \tau_2} \eta_2 \tau_2 \end{pmatrix}, \\ \mathbb{W}(\gamma_1) = \bar{\Gamma}_1, \quad \mathbb{W}(\gamma_2) = \bar{\Gamma}_2, \text{ and } \mathrm{Cov}(\gamma_1, \gamma_2) = \bar{\Gamma}_{12}. \quad (2.14) $$

Then,

$$ \Gamma_1 = \mathrm{diag}(\bar{\Gamma}_1, \dots, \bar{\Gamma}_1), \quad \Gamma_2 = \mathrm{diag}(\bar{\Gamma}_2, \dots, \bar{\Gamma}_2), \text{ and } \Gamma_{12} = \mathrm{diag}(\bar{\Gamma}_{12}, \dots, \bar{\Gamma}_{12}). \quad (2.15) $$
### 2.2.3 EM estimation

Let $\theta$ be the vector of unknown parameters in $\beta_1, \beta_2, \Gamma, \Sigma_1, \Sigma_2$. The EM algorithm requires an initial value of $\theta$ and some expressions (estimators) to update until convergence. In the next two subsections we provide these estimators, their initial values and the stopping criterion.

#### EM-based estimators of parameters

**Theorem 1.** Suppose that $\mathbf{y} = (y_1^\top, y_2^\top)^\top$ satisfies the model based on Eq. (2.2, 2.3 and 2.4) and $\theta$ the vector of its unknown parameters while $\theta_{old}$ is the previous value of $\theta$ provided by the EM algorithm. Let $f(\mathbf{y}, \gamma|\theta)$ be the joint density function of $\mathbf{y}$ and $\gamma$ given $\theta$, and $Q(\theta|\theta_{old}) = \mathbb{E}[\log f(\mathbf{y}, \gamma|\theta)|\mathbf{y}, \theta_{old}]$. Let $M$ be the mapping $\theta_{old} \mapsto M(\theta_{old}) = \hat{\theta}$ such that:

$$M(\theta_{old}) = \arg\max_{\theta} Q(\theta|\theta_{old}) \qquad (2.16)$$

Then, the EM-based estimator of $\theta$, i.e. $\hat{\theta}$, is expressed through:
for $k \in \{1, \dots, m\}$,

$$\hat{\beta}_k = (X_k^\top \Sigma_k^{-1} X_k)^{-1} X_k^\top \Sigma_k^{-1} (y_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \theta_{old}]), \qquad (2.17)$$

$$\hat{\Gamma} = \mathbb{V}[\gamma|\mathbf{y}, \theta_{old}] + \mathbb{E}[\gamma|\mathbf{y}, \theta_{old}] \mathbb{E}[\gamma|\mathbf{y}, \theta_{old}]^\top, \qquad (2.18)$$

$$\hat{\Sigma}_k = Z_k \mathbb{V}[\gamma_k | \mathbf{y}, \theta_{old}] Z_k^\top + (y_k - X_k \beta_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \theta_{old}]) (y_k - X_k \beta_k - Z_k \mathbb{E}[\gamma_k | \mathbf{y}, \theta_{old}])^\top, \quad (2.19)$$

where,

$$\mathbb{E}[\gamma_k|\mathbf{y}, \theta_{old}] = \mathrm{Cov}(\gamma_k, \mathbf{y}|\theta_{old}) \mathbb{V}[\mathbf{y}|\theta_{old}]^{-1} (\mathbf{y} - \mathbb{E}[\mathbf{y}|\theta_{old}]), \qquad (2.20)$$

$$\mathbb{V}[\gamma_k|\mathbf{y}, \theta_{old}] = \mathbf{\Gamma}_k - \mathrm{Cov}(\gamma_k, \mathbf{y}|\theta_{old}) \mathbb{V}[\mathbf{y}|\theta_{old}]^{-1} \mathrm{Cov}(\gamma_k, \mathbf{y}|\theta_{old})^\top \qquad (2.21)$$

and

$$\mathbb{V}[\mathbf{y}|\theta_{old}] = \begin{pmatrix} Z_1 \Gamma_1 Z_1^\top + \Sigma_1 & Z_1 \Gamma_{12} Z_2^\top \\ Z_2 \Gamma_{12}^\top Z_1^\top & Z_2 \Gamma_2 Z_2^\top + \Sigma_2 \end{pmatrix}, \quad \mathrm{Cov}(\gamma, \mathbf{y}|\theta_{old}) = \begin{pmatrix} Z_1 \Gamma_1 & Z_1 \Gamma_{12} \\ Z_2 \Gamma_{12}^\top & Z_2 \Gamma_2 \end{pmatrix}^\top, \qquad (2.22)$$

$$\mathrm{Cov}(\gamma_1, \mathbf{y}|\theta_{old}) = \begin{pmatrix} Z_1 \Gamma_1 \\ Z_2 \Gamma_{12}^\top \end{pmatrix}^\top, \quad \mathrm{Cov}(\gamma_2, \mathbf{y}|\theta_{old}) = \begin{pmatrix} Z_1 \Gamma_{12} \\ Z_2 \Gamma_2 \end{pmatrix}^\top. \qquad (2.23)$$
Proof. For $k \in \{1, \dots, m\}$, $\hat{\beta}_k$, $\hat{\Sigma}_k$ and $\hat{\Gamma}$ optimize the quantity:

$$Q(\theta|\theta_{\text{old}}) = \mathbb{E}[\log f(\mathbf{y}, \gamma|\theta)|\mathbf{y}, \theta_{\text{old}}] \quad (2.24)$$

where $f(\mathbf{y}, \gamma|\theta)$ is the joint density function of the observed data $\mathbf{y}$ and the random effect $\gamma$. In the case of $m=2$, we have:

$$ \begin{aligned} f(\mathbf{y}, \gamma|\theta) &= f(\mathbf{y}|\gamma, \theta)f(\gamma|\theta) \\ &= f(y_1|\gamma_1, \theta)f(y_2|\gamma_2, \theta)f(\gamma|\theta) \\ &= (2\pi)^{-(N_1+N_2+q)/2}|\Sigma_1|^{-1/2}|\Sigma_2|^{-1/2}|\Gamma|^{-1/2} \exp \left\{ -\frac{1}{2} \gamma^\top \Gamma^{-1} \gamma \right. \\ & \qquad \left. -\frac{1}{2}(y_1 - X_1\beta_1 - Z_1\gamma_1)^\top \Sigma_1^{-1} (y_1 - X_1\beta_1 - Z_1\gamma_1) \right. \\ & \qquad \left. -\frac{1}{2}(y_2 - X_2\beta_2 - Z_2\gamma_2)^\top \Sigma_2^{-1} (y_2 - X_2\beta_2 - Z_2\gamma_2) \right\} \end{aligned} \quad (2.25) $$

Since $f$ is a multivariate Gaussian, using the dominated convergence theorem and the derivative under the integral sign, the differential of $Q(\theta|\theta_{\text{old}})$ yields:

$$ \begin{aligned} dQ(\theta|\theta_{\text{old}}) &= \mathbb{E} \left[ -\frac{1}{2} \operatorname{tr} \left\{ \Sigma_1^{-1} d\Sigma_2 + \Sigma_2^{-1} d\Sigma_2 + \Gamma^{-1} d\Gamma \right\} \right. \\ & \qquad \left. -\frac{1}{2} \operatorname{tr} \left\{ -2(y_1 - X_1\beta_1 - Z_1\gamma_1)^\top \Sigma_1^{-1} X_1 d\beta_1 \right. \right. \\ & \qquad \left. -\Sigma_1^{-1} (y_1 - X_1\beta_1 - Z_1\gamma_1)(y_1 - X_1\beta_1 - Z_1\gamma_1)^\top \Sigma_1^{-1} d\Sigma_1 \right\} \\ & \qquad \left. -\frac{1}{2} \operatorname{tr} \left\{ -2(y_2 - X_2\beta_2 - Z_2\gamma_2)^\top \Sigma_2^{-1} X_2 d\beta_2 \right. \right. \\ & \qquad \left. -\Sigma_2^{-1} (y_2 - X_2\beta_2 - Z_2\gamma_2)(y_2 - X_2\beta_2 - Z_2\gamma_2)^\top \Sigma_2^{-1} d\Sigma_2 \right\} \\ & \qquad \left. -\frac{1}{2} \operatorname{tr} \left\{ -\gamma^\top \Gamma^{-1} d\Gamma \Gamma^{-1} \gamma \right\} \bigg|_{\mathbf{y}, \theta_{\text{old}}} \right] \end{aligned} \quad (2.26) $$

$$ \begin{aligned} &= -\frac{1}{2} \operatorname{tr} \left\{ \Sigma_1^{-1} d\Sigma_2 + \Sigma_2^{-1} d\Sigma_2 + \Gamma^{-1} d\Gamma \right\} \\ & \quad + \operatorname{tr}(y_1 - X_1\beta_1 - Z_1\mathbb{E}[\gamma_1|\mathbf{y}, \theta_{\text{old}}])^\top \Sigma_1^{-1} X_1 d\beta_1 \\ & \quad + \frac{1}{2} \operatorname{tr} \left\{ \Sigma_1^{-1} \mathbb{E}[(y_1 - X_1\beta_1 - Z_1\gamma_1)^\top (y_1 - X_1\beta_1 - Z_1\gamma_1)|\mathbf{y}, \theta_{\text{old}}] \Sigma_1^{-1} d\Sigma_1 \right\} \\ & \quad + \operatorname{tr}(y_2 - X_2\beta_2 - Z_2\mathbb{E}[\gamma_2|\mathbf{y}, \theta_{\text{old}}])^\top \Sigma_2^{-1} X_2 d\beta_2 \\ & \quad + \frac{1}{2} \operatorname{tr} \left\{ \Sigma_2^{-1} \mathbb{E}[(y_2 - X_2\beta_2 - Z_2\gamma_2)^\top (y_2 - X_2\beta_2 - Z_2\gamma_2)|\mathbf{y}, \theta_{\text{old}}] \Sigma_2^{-1} d\Sigma_2 \right\} \\ & \quad + \operatorname{tr} \frac{1}{2} \left\{ \Gamma^{-1} \mathbb{E}[\gamma \gamma^\top | \mathbf{y}, \theta_{\text{old}}] \Gamma^{-1} d\Gamma \right\} \end{aligned} \quad (2.27) $$

Partial derivatives of $Q(\theta|\theta_{\text{old}})$ yield:
for $k \in \{1, \dots, m\}$,

$$ \frac{\partial Q(\theta|\theta_{\text{old}})}{\partial \beta_k} = (y_k - X_k\beta_k - Z_k\mathbb{E}[\gamma_k|\mathbf{y}, \theta_{\text{old}}])^\top \Sigma_k^{-1} X_k, $$
$$ \frac{\partial Q(\theta|\theta_{\text{old}})}{\partial \Sigma_k} = -\frac{1}{2} \Sigma_k^{-1} + \Sigma_k^{-1} \mathbb{E}[(y_k - X_k\beta_k - Z_k\gamma_k)^{\top} (y_k - X_k\beta_k - Z_k\gamma_k) | \mathbf{y}, \theta_{\text{old}}] \Sigma_k^{-1} $$

and

$$ \frac{\partial Q(\theta|\theta_{\text{old}})}{\partial \Gamma} = \frac{1}{2} (\Gamma^{-1} \mathbb{E}[\gamma \gamma^{\top} | \mathbf{y}, \theta_{\text{old}}] \Gamma^{-1} - \Gamma^{-1}). $$

We then get EM-based estimators by setting these partial derivatives equal to zero. $\mathbb{E}[\gamma_k|\mathbf{y}, \theta_{\text{old}}]$ and $\mathbb{V}[\gamma_k|\mathbf{y}, \theta_{\text{old}}]$ are straightforward to get since $(\gamma_k^{\top}, \mathbf{y}^{\top})^{\top}$ is a multivariate Gaussian. $\square$

**Initialization and stopping criterion of the algorithm** Various ways exist for obtaining starting values for $\hat{\beta}_k$, $\hat{\Gamma}$, $\hat{\Sigma}_k$, for $k = 1, \dots, m$. Taking inspiration from [103] and [172], we have separately fitted each dimension of the model by using the lme4 package [15] under the R software and have used marginal estimated parameters to initialize $\hat{\beta}_k$ and $\hat{\Sigma}_k$. We then keep the expected random effects $\tilde{\gamma} = \hat{\mathbb{E}}[\gamma|\mathbf{y}]$ to initialize $\hat{\Gamma}$ by

$$ \frac{1}{n-1} \sum_{i=1}^{n} \tilde{\gamma}_i \tilde{\gamma}_i^{\top} \qquad (2.28) $$

The stopping criterion is related to the relative error of the components of $\theta$ as follows:

$$ \max_j \left| \frac{\theta_j^{(r)} - \theta_j^{(r+1)}}{\theta_j^{(r+1)}} \right| < \text{tol} \qquad (2.29) $$

where $(r)$ is the iteration index and $\theta_j$ the $j$th component of $\theta$. $\text{tol} = 10^{-5}$ seems to work well in practice.

### 2.2.4 Test of the significance of $\hat{\text{Cor}}(\gamma_1, \gamma_2)$

After the calculation of $\Gamma$ on dataset, we sometimes need to investigate if the correlation between marginal random effects is statistically significant, by testing $H_0$: $\text{Cor}(\gamma_1, \gamma_2) = 0$ against $H_1$: $\text{Cor}(\gamma_1, \gamma_2) \neq 0$. The result of this test can help to decide if the bivariate analysis is justified or not. We perform the likelihood ratio (LR) test to choose between $H_0$ and $H_1$. We calculated $S$, the statistic of the likelihood ratio test.

$$ S = -2 \log \left( \frac{L(\theta|\text{data}, H_0)}{L(\theta|\text{data}, H_1)} \right) \qquad (2.30) $$

where $L(\theta|\text{data}, H_0)$ and $L(\theta|\text{data}, H_1)$ are the likelihood of $\theta$ under $H_0$ and $H_1$, respectively. Under suitable and standard conditions, $S \sim \chi^2(df)$, asymptotically, under $H_0$ [200]. With $df$ the difference in the number of parameters between $L(\theta|\text{data}, H_0)$ and $L(\theta|\text{data}, H_1)$.
## 2.3 Results and Discussion

### 2.3.1 Simulation studies

In this section, simulation studies are used to investigate the computational properties of the EM-based estimators. For the sake of simplicity, these simulation studies are conducted using simulated bivariate longitudinal data sets. Through these studies, we pursue two objectives: the first is to assess the accuracy of parameter estimates and the second is to analyze the power of the likelihood ratio test performed via these EM-based estimators. In the following paragraph, we explain how we choose the parameters that have been used to simulate the working longitudinal data sets.

**The working data sets** We suppose that we are following up a sample of subjects where the goal is to evaluate how the growth of the weight and the height of the individuals of this population are jointly explained by the sex, the score of nutrition (Nscore) and the age. We randomly choose through a uniform distribution the score of nutrition between 20 and 50, and the age between 18 and 37, using the R software. All the analysis in this paper are done using the R software. The subject's sex is also randomly chosen. The model under which we simulate the data sets is the following:

$n$ indicating the total number of subjects, for $i = 1, \dots, n$

$$ 
\begin{aligned}
\text{weight}_i &= (\mathbb{1}_{n_i}, \text{sex}_i, \text{Nscore}_i, \text{age}_i)\beta_1 + (\mathbb{1}_{n_i}, \text{Nscore}_i)\gamma_{1i} + \varepsilon_{1i} \\
\text{height}_i &= (\mathbb{1}_{n_i}, \text{sex}_i, \text{Nscore}_i, \text{age}_i)\beta_2 + (\mathbb{1}_{n_i}, \text{Nscore}_i)\gamma_{2i} + \varepsilon_{2i}
\end{aligned}
\quad (2.31) $$

with

$$ \gamma_i = \begin{pmatrix} \gamma_{1i} \\ \gamma_{2i} \end{pmatrix} \sim \mathcal{N}(0, \bar{\Gamma}), \quad \varepsilon_{1i} \sim \mathcal{N}(0, \sigma_1^2 \mathbf{I}_{n_i}), \quad \varepsilon_{2i} \sim \mathcal{N}(0, \sigma_2^2 \mathbf{I}_{n_i}), \quad \gamma_i \perp \varepsilon_{1i} \perp \varepsilon_{2i} \quad (2.32) $$

The random effect related to the dependent variable 'weight' or 'height' is a vector composed by one random intercept and one random slope in the direction of the covariate 'Nscore'. The total number of observations is denoted by $N$.

We randomly choose $\beta_1, \beta_2, \sigma_1$ and $\sigma_2$ whose values are in the first column of Table 2.2. $\bar{\Gamma}$ is also randomly chosen such that it is positive definite, with the following form:

$$ \bar{\Gamma} = \begin{pmatrix} \eta_1^2 & \rho\eta_1\eta_2 & \rho\eta_1\tau_1 & \rho\eta_1\tau_2 \\ \rho\eta_1\eta_2 & \eta_2^2 & \rho\eta_2\tau_1 & \rho\eta_2\tau_2 \\ \rho\eta_1\tau_1 & \rho\eta_2\tau_1 & \tau_1^2 & \rho\tau_1\tau_2 \\ \rho\eta_1\tau_2 & \rho\eta_2\tau_2 & \rho\tau_1\tau_2 & \tau_2^2 \end{pmatrix} \quad (2.33) $$

The covariance between the random effects $\gamma_1$ and $\gamma_2$ is set, intentionally,

$$ \text{Cov}(\gamma_1, \gamma_2) = \rho \begin{pmatrix} \eta_1 \tau_1 & \eta_1 \tau_2 \\ \eta_2 \tau_1 & \eta_2 \tau_2 \end{pmatrix} \quad (2.34) $$
<table><caption>Table 2.2: Comparative table of true values of parameters and estimates based on 1000 repli-cations using true values of parameters.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Empirical mean</th><th>Empirical Sd</th><th>Bias</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50.67</td><td>50.669</td><td>0.763</td><td>0.000</td></tr><tr><td>-4.80</td><td>-4.779</td><td>0.811</td><td>0.021</td></tr><tr><td>14.00</td><td>14.012</td><td>0.345</td><td>0.012</td></tr><tr><td>2.70</td><td>2.700</td><td>0.016</td><td>0.000</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>13.20</td><td>13.263</td><td>1.077</td><td>0.063</td></tr><tr><td>-2.80</td><td>-2.796</td><td>1.186</td><td>0.003</td></tr><tr><td>27.00</td><td>27.000</td><td>0.068</td><td>0.000</td></tr><tr><td>1.70</td><td>1.699</td><td>0.019</td><td>0.000</td></tr><tr><td>σ<sub>1</sub></td><td>5.80</td><td>5.796</td><td>0.062</td><td>0.003</td></tr><tr><td>σ<sub>2</sub></td><td>7.60</td><td>7.602</td><td>0.082</td><td>0.002</td></tr></tbody></table>

in order to be able to decrease or increase the correlation between the marginal random effects $\gamma_1$ and $\gamma_2$, by changing the value of $\rho$, without losing the positive definiteness of $\bar{\Gamma}$. This property of $\bar{\Gamma}$ will be used to assess the power of the likelihood ratio test through simulations, by changing the value of $\rho$. We simulate 1000 data sets with $\rho = 0.8$, in order to assess the accuracy of estimates using the EM-based estimators. With $\rho = 0.8$, the randomly chosen $\bar{\Gamma}$ is

$$ \bar{\Gamma} = \begin{pmatrix} 27.77 & 18.80 & 41.70 & 4.93 \\ 18.80 & 36.00 & 47.47 & 5.62 \\ 41.70 & 47.47 & 97.81 & 8.91 \\ 4.93 & 5.62 & 8.91 & 1.37 \end{pmatrix} \qquad (2.35) $$

**Empirical accuracy of the estimates** The 1000 data sets simulated in order to assess the accu- racy of the estimates performed using the EM-based estimators contain $N = 5000$ observations provided by $n = 300$ independent subjects.

The mean and the standard deviation of the 1000 estimates are presented, respectively, in the second and the third column of the Table 2.2. The bias of the parameter estimates, which is the absolute difference between the true value of the parameter and the mean of the 1000 estimates, is calculated as measure of performance. These bias are contained in the forth column of the Table 2.2.

$\hat{\mathbb{E}}[\hat{\bar{\Gamma}}]$, $\hat{\sigma}_{\hat{\bar{\Gamma}}}$ and $\text{Bias}(\hat{\bar{\Gamma}})$ (Equations 2.36, 2.37) contain, respectively, the empirical mean, the empirical standard deviation and the empirical bias of $\bar{\Gamma}$.

$$ \hat{\mathbb{E}}[\hat{\bar{\Gamma}}] = \begin{pmatrix} 27.86 & 18.71 & 41.19 & 4.93 \\ 18.71 & 35.73 & 47.09 & 5.58 \\ 41.19 & 47.09 & 95.93 & 8.87 \\ 4.93 & 5.58 & 8.87 & 1.36 \end{pmatrix}; \quad \hat{\sigma}_{\hat{\bar{\Gamma}}} = \begin{pmatrix} 5.32 & 2.99 & 6.37 & 0.63 \\ 2.98 & 2.86 & 5.11 & 0.50 \\ 6.37 & 5.11 & 13.73 & 0.97 \\ 0.63 & 0.50 & 0.97 & 0.11 \end{pmatrix} \qquad (2.36) $$
<figure><img src="image_3.png" /><figcaption>Figure 2.1: Coefficients of variation of entries of Γ. N = 5000 observations and n = 300 subjects.</figcaption></figure>

$$ \text{Bias}(\hat{\Gamma}) = \begin{pmatrix} 0.085 & 0.090 & 0.515 & 0.001 \\ 0.089 & 0.265 & 0.382 & 0.045 \\ 0.514 & 0.382 & 1.881 & 0.039 \\ 0.001 & 0.045 & 0.039 & 0.007 \end{pmatrix} \qquad (2.37) $$

The bias contained in the estimates of $\beta_k$ and $\sigma_k$ ranges from 0.000 to 0.063 (Table 2.2), and the bias contained in the estimates of the entries of $\bar{\Gamma}$ ranges from 0.001 to 1.881 (Equation 2.37). These results show that $\hat{\beta}_k$ and $\hat{\sigma}_k$ (i.e. $\hat{\Sigma}_k$) seem unbiased when $\hat{\Gamma}$ is biased.

The estimates of $\bar{\Gamma}$ appear to be poorer than the estimates of all other parameters. In order to investigate which entries of $\bar{\Gamma}$ are particularly poorly estimated, we calculate the coefficients of variation (CV) of these entries. The CV computed here is obtained by dividing the standard deviation of the estimates by the true value of each entry of $\bar{\Gamma}$. The CVs give an idea of the variability of estimates around the true values and enable to compare these variabilities between them. A particularly large value of CV could lead us to suspect that the corresponding input is particularly poorly estimated. Here, the CV ranges from 0.08 to 0.19, and is represented by the Fig 2.1 for more visibility. Given these CV values, it seems that none of the entries of $\bar{\Gamma}$ is particularly poorly estimated.

**Deep investigation on the estimates' accuracy** Here, we compute the Mean Square Error (MSE) of the EM-based estimators with N = 600, 1000 and N = 3000 across n = 50, 60, 100 and 300 to investigate how both values of n and N affect the quality of the estimates. For each value of n
<table><caption>Table 2.3: Mean Square Error of EM-based estimator with 95% CI estimated on 1000 repli-cations for various values of *n* and *N*.</caption><thead><tr><th>Parameter</th><th>n</th><th>N = 600</th><th>N = 1000</th><th>N = 3000</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50</td><td>1.27 (0.01 - 4.64)</td><td>1.10 (0.00 - 3.83)</td><td>0.71 (0.00 - 2.57)</td></tr><tr><td>60</td><td>1.38 (0.00 - 5.56)</td><td>0.99 (0.00 - 3.90)</td><td>0.65 (0.00 - 2.44)</td></tr><tr><td>100</td><td>1.31 (0.00 - 4.79)</td><td>0.87 (0.00 - 3.41)</td><td>0.46 (0.00 - 1.58)</td></tr><tr><td>300</td><td>1.64 (0.00 - 5.93)</td><td>0.80 (0.00 - 3.23)</td><td>0.29 (0.00 - 1.13)</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>50</td><td>2.13 (0.01 - 7.22)</td><td>1.73 (0.00 - 6.21)</td><td>1.08 (0.00 - 4.17)</td></tr><tr><td>60</td><td>2.12 (0.00 - 8.06)</td><td>1.63 (0.01 - 6.66)</td><td>0.95 (0.00 - 3.52)</td></tr><tr><td>100</td><td>1.88 (0.00 - 7.14)</td><td>1.29 (0.00 - 4.76)</td><td>0.62 (0.00 - 2.24)</td></tr><tr><td>300</td><td>2.30 (0.00 - 8.54)</td><td>1.23 (0.00 - 4.57)</td><td>0.43 (0.00 - 1.65)</td></tr><tr><td rowspan="4">σ<sub>1</sub></td><td>50</td><td>0.03 (0.00 - 0.10)</td><td>0.10 (0.00 - 0.07)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.07)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>100</td><td>0.04 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.08)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td>300</td><td>0.06 (0.00 - 0.22)</td><td>0.03 (0.00 - 0.11)</td><td>0.01 (0.00 - 0.02)</td></tr><tr><td rowspan="4">σ<sub>2</sub></td><td>50</td><td>0.05 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>60</td><td>0.06 (0.00 - 0.21)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>100</td><td>0.06 (0.00 - 0.24)</td><td>0.04 (0.00 - 0.15)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>300</td><td>0.09 (0.00 - 0.36)</td><td>0.04 (0.00 - 0.18)</td><td>0.01 (0.00 - 0.05)</td></tr><tr><td rowspan="4">Γ̄</td><td>50</td><td>400.91 (1.52 - 1274.32)</td><td>670.46 (2.36 - 2536.73)</td><td>489.07 (2.02 - 1840.60)</td></tr><tr><td>60</td><td>706.45 (2.58 - 2497.42)</td><td>620.34 (2.70 - 2293.58)</td><td>408.79 (0.85 - 1597.77)</td></tr><tr><td>100</td><td>701.50 (3.79 - 2747.70)</td><td>477.55 (2.08 - 1839.26)</td><td>283.25 (1.97 - 1023.77)</td></tr><tr><td>300</td><td>798.28 (3.65 - 2603.65)</td><td>547.83 (3.03 - 1721.46)</td><td>199.54 (0.81 - 736.08)</td></tr></tbody></table>

and N, we simulate 1000 data sets on which we estimate the model parameters and compute the MSE of these estimates.

Without surprise, Table 2.3 shows that the quality of estimates is clearly improved when both *n* and *N* grow. Estimations performed on dataset containing *N* = 3000 observations are more accurate than those performed with *N* = 600, observing the maximum value of the MSE in each case. For *N* = 600, information contained in Table 2.3 shows that the MSE related to *n* = 60 (60 subjects) are better than those related to *n* = 300. This result suggests a good tradeoff between the number of subjects and the total number of observations in order to have accurate estimates, especially if the number of observations is not very high. Once again, it appears that $\hat{\Gamma}$ (Table 2.3) has the highest MSE for all values of *n* and *N*.

**The bivariate likelihood ratio test** Considering the random effects covariance matrix $\bar{\Gamma}$ (see Equation 3.44), the related correlation matrix is

$$ \begin{pmatrix} 1.00 & 0.59 & 0.80 & 0.80 \\ 0.59 & 1.00 & 0.80 & 0.80 \\ 0.80 & 0.80 & 1.00 & 0.77 \\ 0.80 & 0.80 & 0.77 & 1.00 \end{pmatrix}. \qquad (2.38) $$

That is, the matrix of the correlations between the marginal random effects (i.e., the random effects related to the two dependent variables) is

$$ \operatorname{Cor}(\gamma_1, \gamma_2) = \begin{pmatrix} 0.80 & 0.80 \\ 0.80 & 0.80 \end{pmatrix} \qquad (2.39) $$
whereas the estimate (on one of the previous simulated data) of this matrix, Cor(γ₁, γ₂), of the correlations between the marginal random effects, is

$$ \widehat{\text{Cor}}(\gamma_1, \gamma_2) = \begin{pmatrix} 0.77 & 0.78 \\ 0.90 & 0.74 \end{pmatrix} \qquad (2.40) $$

If we decide to test $H_0$: Cor(γ₁, γ₂) = 0 against $H_1$: Cor(γ₁, γ₂) ≠ 0 in the case of these simulated data, we must know the distribution of the LR statistic S. In order to approximate the distribution of S, under $H_0$, we proceed to an extensive simulation study in the next paragraph.

**Empirical distribution of S under $H_0$** In this paragraph, our goal is to investigate about the empirical law of the LR statistic S, under $H_0$, when the size N of the data set increases. The simulated data sets used in this paragraph are also of bivariate longitudinal type, with N the total number of observations coming from n subjects. We choose N as an arithmetic sequence ranging from 50 to 2000, where the common difference is 50. We choose n = N/5 as it is sufficient to have two observations per subject for fitting the model. When N/n = 1, the random-effects parameters and the residual variance are unidentifiable [140].

The expected (standard) asymptotic distribution of S, under $H_0$, is a $\chi^2(4)$. This may be explained by the fact that Cov(γ₁, γ₂) and its transpose, Cov(γ₁, γ₂)ᵀ, contain four entries, respectively, and $\bar{\Gamma}$ contains Cov(γ₁, γ₂) and Cov(γ₁, γ₂)ᵀ. Therefore, the difference between the number of entries of $\bar{\Gamma}$ which need to be estimated with $L(\theta|$data, $H_0$) and $L(\theta|$data, $H_1$), respectively, is df = 4. Precisely, the parameters of interest are $\rho_{\eta_1\tau_1}$, $\rho_{\eta_1\tau_2}$, $\rho_{\eta_2\tau_1}$ and $\rho_{\eta_2\tau_2}$ (see Equation 2.14).

<figure><img src="image_4.png" /><figcaption>Figure 2.2: Empirical analysis of the asymptotic distribution of the LR statistic S under $H_0$, using longitudinal data sets (200 replications) with size $N \in \{50, 100, 150, 200, \dots, 2000\}$ coming from $n \in \{10, 20, 30, \dots, 400\}$ subjects. An asymptotic distribution of $\chi^2(4)$ is assumed and the Kolmogorov-Smirnov test's p-value (at log10 scale) is plotted against the total number of observations of the data set that has served to compute the LR statistic S. The blue curve is obtained by applying the empirical Bartlett correction to S and the red curve is obtained without correction. The horizontal dashed line represents $log10(0.05)$.</figcaption></figure>
Fig 2.2 assumes an asymptotic distribution of $\chi^2(4)$ and plots the Kolmogorov-Smirnov test's p-value (at log10 scale) against the total number of observations of the data set that has served to compute the LR statistic $S$. The blue curve is obtained by applying the empirical Bartlett correction to $S$ and the red curve is obtained without correction. The horizontal dashed line represents $\log10(0.05)$. The empirical Bartlett corrected $S$, say $\hat{S}_B$, can be expressed as $\hat{S}_B = df \times S/\hat{E}[S|H_0]$. This Bartlett correction is applied in order to avoid the small size distortion of the $\chi^2(df)$ distribution, when performing the LR test using a data set of small size [12]. Fig 2.2 thus helps to investigate how the LR distribution performs in finite and small dimension. It also helps to investigate, in the case of this bivariate correlation test, how the Bartlett correction helps to avoid the small size distortion of the chi-square approximation. As the total number of observations increases, the curves (red and blue) reach the dashed line, gradually. Assuming the $\chi^2(4)$ distribution of $S$, it seems important to work with a data set containing at least 500 observations coming from at least 2 subjects, and to apply the Bartlett correction in order to avoid the breakdown of the procedure.

<figure><img src="image_5.png" /><figcaption>Figure 2.3: Empirical analysis of the asymptotic distribution of the LR statistic $S$ under $H_0$, using longitudinal data sets (200 replications) with size $N \in \{50, 100, 150, 200, \dots, 2000\}$ coming from $n \in \{10, 20, 30, \dots, 400\}$ subjects. An asymptotic distribution of $\chi^2(4)$ is assumed and the type I error (at log10 scale) is plotted against the total number of observations of the data set that has served to compute the LR statistic $S$. The blue curve is obtained by applying the empirical Bartlett correction to $S$ and the red curve is obtained without correction. The horizontal dashed lines represent the significance levels of 5% and 10%, respectively.</figcaption></figure>

The type I error is generally controlled by the significance level of 10% (red and blue curves of Fig 2.3). It is clear that the control is almost full with the Bartlett correction (blue curve of Fig 2.3).

By simulating 1000 × 3000 realizations of $\chi^2(4)$ distribution, we plot the red sheath represented in Fig 2.4. This sheath corresponds to the minimum and the maximum of the simulated $\chi^2(4)$ realizations. The blue curve inside the red sheath represents the empirical LR statistics obtained from the 3000 simulated data sets under $H_0$. This figure (Fig 2.4) shows that the asymptotic
<figure><img src="image_6.png" /><figcaption>Figure 2.4: Empirical analysis of the asymptotic distribution of the LR statistic S under $H_0$, using 3000 (replications) simulated longitudinal data sets (under $H_0$) of size $N = 15000$ coming from $n = 500$ subjects. The minimum and the maximum of 1000 × 3000 simulated realizations of $\chi^2(4)$ are used to construct the red sheath. The blue curve represents the LR statistics related to the bivariate correlation test.</figcaption></figure>

distribution of LR statistic related to the bivariate correlation test is not violated, since the blue curve does not go out of the red sheath.

**Empirical power of the bivariate correlation test** In order to analyze the power of this likelihood ratio test performed with EM-based estimates, we calculate S on data sets which have been simulated under $H_0$ and $H_1$, respectively, leading to what we named S0 and S1 vectors containing the resulting values of S. We then plot a ROC curve with S0 and S1, where S0 is the vector of the cases while S1 contains the controls. We calculate S0 and S1 in different situations where we have changed the value of $\rho$ in the following configuration:

$$Cov(\gamma_1, \gamma_2) = \rho \begin{pmatrix} \eta_1 \tau_1 & \eta_1 \tau_2 \\ \eta_2 \tau_1 & \eta_2 \tau_2 \end{pmatrix} \qquad (2.41)$$

We maintain fixed $\eta_1 = 5.27$, $\eta_2 = 6.00$, $\tau_1 = 9.89$, $\tau_2 = 1.17$ and change $\rho$ ($\in \{0.1, 0.2, 0.3, \dots, 0.9\}$). The number of subjects ($n$) and the total number of observations ($N$) have also been modified throughout these simulation studies. In each case, the estimated Area Under Curve (AUC) of the ROC curve with its confidence interval have been recorded to produce Fig 2.5.

With $n = 50$ subjects, we detect, indeed, a correlation of 0.6 when the total number of observations is $N = 3000$; in contrast, if the total number of observations is $N = 600$, we perfectly detect a correlation of 0.7.

Unsurprisingly, confidence intervals of AUC are also more accurate with $N = 3000$ than with $N = 600$. With a sufficient number of observations and subjects, weak correlations are easily detected. For example, we perfectly detect a correlation of 0.2 with $N = 3000$ and $n = 300$ where
<figure><img src="image_7.png" /><figcaption>Figure 2.5: **Empirical analysis of the power of the correlation test.** AUC values of ROC curves with their confidence interval computed for different $\rho$, number of subjects ($n$) and observations ($N$). Left panel for $N = 600, n = 50, 60, 100$. Right panel for $N = 3000, n = 50, 100, 300$.</figcaption></figure>

AUC = 0.96(0.93 – 0.98) according to Fig 2.5. However, we detect quite well a correlation of 0.3 with $N = 600$ and $n = 60$ where AUC = 0.81(0.73 – 0.88).

In the case where estimates are of a higher quality (because they are performed on data sets having a sufficient number of observations $N = 5000$ and subjects $n = 300$), we plot ROC curves with low values of $\rho$ (0.1, 0.2 and 0.3). We then show in Fig 2.6, the estimated AUC and its 95% confidence interval.

Fig 2.6 shows that EM-based estimators lead to a good power of the bivariate correlation test, when we have a sufficient number of observations and subjects in the longitudinal study case. This goodness of the power of the bivariate correlation test persists when the correlation between marginal random effects is low (about 0.2).

### 2.3.2 Applications on real data sets

In this section we analyze two data sets by using the likelihood ratio test through the EM-based estimators presented above. The first dataset is of multivariate multilevel type and the second is, specifically, of longitudinal multivariate type.

**Application to education data in the Netherlands** The data used here are named 'bdf' under the package nlme[142] of the R software. These data contain $N = 3776$ Grade eight students (aged about eleven years) in $n = 208$ elementary schools in the Netherlands [22]. These pupils were tested twice (with an interval of one year between grade seven and grade eight) on their proficiency in Dutch language and arithmetics, where the goal was to investigate which characteristics of schools can account for the differences in the effectiveness of schools with regard to pupils' progress in language and arithmetics. Most of the previous analyses of this dataset were concerned with investigating how the language test score depends on the pupil's intelligence, his
<figure><img src="image_8.png" /><figcaption>Figure 2.6: Analysis of the power of the likelihood ratio test performed under EM-estimators. ROC curves with $\rho \in \{0.1, 0.2, 0.3\}$. N = 5000 observations, n = 300 subjects, 95% CI on AUC.</figcaption></figure>

family's socio-economic status and on related class or school variables. By fitting two independent (separate) models, [22] found that variables in Table 2.4 have a significant effect on post-test scores (language post-test and arithmetic post-test). These variables are: socio-economic status, intelligence score, age, gender and nationality. They also found a significant random slope related to the language pre-test and to the gender in the language post-test model.

Based on these results from [22] and some of their data (n = 131 schools, N = 2287 pupils; age and ethnicity are not present), we have fitted the bivariate linear mixed-effect model where post-test scores are the response variables and covariates are the pre-test scores, socio-economic status, intelligence score, gender and minority (a factor indicating if the pupil is a member of a

<table><caption>Table 2.4: Modeling of covariates on post-test achievement in language and arithmetic from [22]</caption><thead><tr><th>Model</th><th>Language (post)</th><th>Arithmetic (post)</th></tr></thead><tbody><tr><td>Language pre-test</td><td>0.567</td><td>-</td></tr><tr><td>Arithmetic pre-test</td><td>-</td><td>0.413</td></tr><tr><td>Socio-economic status</td><td>0.143</td><td>0.132</td></tr><tr><td>Intelligence score</td><td>0.124</td><td>0.270</td></tr><tr><td>Age</td><td>-0.069</td><td>-0.081</td></tr><tr><td>Gender (female)</td><td>0.187</td><td>-0.103</td></tr><tr><td>Ethnicity (foreign)</td><td>n.s.</td><td>-0.105</td></tr></tbody></table>
<table><caption>Table 2.5: Model configuration</caption><thead><tr><th rowspan="2">Model parts</th><th colspan="2">Variables</th></tr><tr><th>Language (post)</th><th>Arithmetic (post)</th></tr></thead><tbody><tr><td rowspan="4">Fixed effects</td><td>Language (pre)</td><td>Arithmetic (pre)</td></tr><tr><td>Socio-eco. status</td><td>Socio-eco. status</td></tr><tr><td>Intelligence score</td><td>Intelligence score</td></tr><tr><td>Gender</td><td>Gender</td></tr><tr><td>Minority</td><td>Minority</td></tr><tr><td>Random effects (school level)</td><td>1</td><td>1</td></tr><tr><td></td><td>Language (pre)</td><td>Arithmetic (pre)</td></tr></tbody></table>

<table><caption>Table 2.6: Estimated fixed effects and residual standard deviations in the joint bivariate model fitted to school data.</caption><thead><tr><th rowspan="3">Covariates</th><th colspan="4">Response variables</th></tr><tr><th colspan="2">Language (post)</th><th colspan="2">Arithmetic (post)</th></tr><tr><th>Estimate</th><th>p-value</th><th>Estimate</th><th>p-value</th></tr></thead><tbody><tr><td>Intercept</td><td>4.690</td><td></td><td>-1.456</td><td></td></tr><tr><td>Language (pre)</td><td>0.795</td><td>0.000</td><td>-</td><td>-</td></tr><tr><td>Arithmetic (pre)</td><td>-</td><td>-</td><td>0.807</td><td>0.000</td></tr><tr><td>Socio-eco. status</td><td>0.101</td><td>0.000</td><td>0.089</td><td>0.000</td></tr><tr><td>Intelligence score</td><td>0.474</td><td>0.000</td><td>0.810</td><td>0.000</td></tr><tr><td>Gender (female)</td><td>1.778</td><td>0.000</td><td>-0.522</td><td>0.003</td></tr><tr><td>Minority (yes)</td><td>-0.302</td><td>0.589</td><td>-0.545</td><td>0.194</td></tr><tr><td>σ<sub>1</sub> and σ<sub>2</sub></td><td>5.356</td><td>4.066</td><td></td><td></td></tr></tbody></table>

minority group). Random intercepts and random slopes related to pre-test scores are integrated to the model on the school level in the configuration shown by the Table 2.5.

Table 2.6 contains estimated fixed effects and residual standard deviations of the model.

The estimated covariance matrix $\hat{\Gamma}$ of the random effects is:

$$ \hat{\Gamma} = \begin{pmatrix} 15.15 & -0.54 & 18.20 & -0.27 \\ -0.54 & 0.02 & -0.65 & 0.01 \\ 18.20 & -0.65 & 30.69 & -0.48 \\ -0.27 & 0.01 & -0.48 & 0.01 \end{pmatrix} $$

The null hypothesis, $H_0$, that the arithmetic post-test score and the language post-test score are independent, is rejected with a p-value of $1.436 \times 10^{-7}$. This result justifies a joint analysis of
post-scores conditionally on the covariates present in the model and is therefore a supplementary information obtained from the data. The estimated correlation matrix of the random effects is:

$$ \hat{\rho} = \begin{pmatrix} 1.00 & -1.00 & 0.84 & -0.82 \\ -1.00 & 1.00 & -0.84 & 0.82 \\ 0.84 & -0.84 & 1.00 & -1.00 \\ -0.82 & 0.82 & -1.00 & 1.00 \end{pmatrix} $$

Table 2.6 shows that covariates which are significant in the independent models are also significant in the joint model. The Minority covariate is neither significant in the joint model, nor significant in the independent models. $\hat{\Gamma}_{ij}$ and $\hat{\rho}_{ij}$ identify the item which is at the intersection of the row $i$ and the column $j$ of the matrices $\hat{\Gamma}$ and $\hat{\rho}$, respectively. These matrices are filled from the top to the bottom in the order of ($\text{Intercept}_{y_1}$, $\text{Slope}_{y_1}$, $\text{Intercept}_{y_2}$, $\text{Slope}_{y_2}$).

There is a clear inter-school variability with respect to the post-test scores ($\hat{\Gamma}_{11} = 15.15$ and $\hat{\Gamma}_{33} = 30.69$). Everything else being equal, schools that have good scores in arithmetics also have good scores in language ($\hat{\rho}_{31} = 0.84$). The schools in which the differential effect of the pre-test score in arithmetics on the post-test score is strongly negative are in average above the average post-test score in language ($\hat{\rho}_{41} = -0.82$); same as with the language pre-test score ($\hat{\rho}_{32} = -0.84$). This confirms that the scores in language and in arithmetics vary in the same direction in schools. The schools in which the differential effect of the pre-test score (in arithmetics or language) on the post-test score is strongly negative are in average above the average post score ($\hat{\rho}_{21} = -1$ and $\hat{\rho}_{43} = -1$), and vice versa. These schools have strived to bring the level of all pupils above the average. In contrast, pupils with a good initial level maintain their level without becoming excellent. The differential effect of the pre-test score has a very weak variability (arithmetics score: $\hat{\Gamma}_{22} = 0.02$; language score: $\hat{\Gamma}_{44} = 0.01$) and this implies that the pre-test score explains about 0.15% (in arithmetics) and 0.03% (in language) of inter-school variability of post-test scores.

We have fitted the bivariate model without random slopes (with random intercept only) to investigate if it fits more to the data than the model with random slopes, due to the weak variability of these random slopes. The results are presented in Table 2.7, where the estimated fixed effects and their significance generally remain the same.

The estimated covariance matrix, related to results contained in Table 2.7 of random effects is

$$ \hat{\Gamma} = \begin{pmatrix} 5.15 & 5.41 \\ 5.41 & 7.07 \end{pmatrix} \qquad (2.42) $$

which indicates a correlation of $\hat{\rho} = 0.895$ between the random marginal intercepts, confirming a strong positive correlation between post-test scores in arithmetics and language. With a p-value of $8.505 \times 10^{-5}$, the likelihood ratio test indicates that the data are more compatible with the model incorporating random intercepts and random slopes at a time.

Fixed effects seem very strong and do not significantly change between the independent and bivariate models. In contrast, a posterior distribution of random effects changes significantly between the independent model and the joint bivariate model. For example, we plot the joint distribution of random effects conditional on the data concerning School 47 in the education dataset under the
<table><caption>Table 2.7: Results of the model with random intercepts only.</caption><thead><tr><th rowspan="3">Covariates</th><th colspan="4">Response variables</th></tr><tr><th colspan="2">Language (post)</th><th colspan="2">Arithmetic (post)</th></tr><tr><th>Estimate</th><th>p-value</th><th>Estimate</th><th>p-value</th></tr></thead><tbody><tr><td>Intercept</td><td>4.698</td><td></td><td>-1.446</td><td></td></tr><tr><td>Language (pre)</td><td>0.789</td><td>0.000</td><td>-</td><td>-</td></tr><tr><td>Arithmetic (pre)</td><td>-</td><td>-</td><td>0.789</td><td>0.000</td></tr><tr><td>Socio-eco. status</td><td>0.103</td><td>0.000</td><td>0.093</td><td>0.000</td></tr><tr><td>Intelligence score</td><td>0.480</td><td>0.000</td><td>0.809</td><td>0.000</td></tr><tr><td>Gender (female)</td><td>1.788</td><td>0.000</td><td>-0.526</td><td>0.002</td></tr><tr><td>Minority (yes)</td><td>-0.391</td><td>0.489</td><td>-0.498</td><td>0.249</td></tr><tr><td>σ<sub>1</sub> and σ<sub>2</sub></td><td>5.384</td><td>4.091</td><td></td><td></td></tr></tbody></table>

independent model and the joint bivariate model. Fig 2.7a shows the joint posterior distribution of random intercepts under the independent model whereas Fig 2.7b presents the same posterior distribution under the joint bivariate model. A clear difference appears between these two distributions. We notice the same difference between distributions of random intercepts and slopes as shown in Fig 2.7c and Fig 2.7d as well as the joint distribution of random slopes in Fig 2.7e and Fig 2.7f. The joint bivariate model seems to fit more to the present data and we retain it for their analysis.

**Application to malaria immune response data in Benin** The data come from a study which was conducted in nine villages (Avamé centre, Gbédjougo, Houngo, Anavié, Dohinoko, Gbétaga, Tori Cada Centre, Zébè and Zoungoudo) of Tori Bossito area (Southern Benin), where *P. falci-parum* is the most common species in the study area (95%) [47] from June 2007 to January 2010. The aim of this study was to evaluate the determinants of malaria incidence in the first months of life of child in Benin. Details of the follow-up procedures have been published elsewhere [108].

**Data description** Mothers (n = 620) were enrolled at delivery and their newborns were actively followed-up during the first year of life. One questionnaire was conducted to gather information on women's characteristics (age, parity, use of Intermittent Preventive Treatment during pregnancy (IPTp) and bed net possession) and on the course of their current pregnancy. Material peripheral blood as well as cord blood were collected into Vacutainer® EDTA (Ethylene diaminetetraacetic acid) tubes. At birth, newborn's weight and length were measured by midwives and gestational age was estimated using the Ballard method [9].

During the follow-up of newborns, axillary temperature was measured weekly. Symptomatic malaria cases, defined as fever (> 37.5°C) with TBS and/or RDT positive, were treated with an artemisinin-based combination therapy as recommended by the Benin National Malaria Control Program. Systematically, TBS were made every month to detect asymptomatic infections. Every three months, venous blood was sampled to quantify the level of antibody against malaria promised candidate vaccine antigens. The environmental risk of exposure to malaria was modeled for each
<table><caption>Table 2.8: Variables present in the analyzed dataset</caption><thead><tr><td>Variable</td><td>Description</td></tr></thead><tbody><tr><td>id</td><td>Child ID</td></tr><tr><td>conc.Y</td><td>concentration of Y</td></tr><tr><td>conc_CO.Y</td><td>Measured concentration of Y in the umbilical cord blood</td></tr><tr><td>conc_M3.Y</td><td>Predicted concentration of Y in the child's peripheral blood at 3 months</td></tr><tr><td>ap</td><td>Placental apposition</td></tr><tr><td>hb</td><td>Hemoglobin level</td></tr><tr><td>inf_trim</td><td>Number of malaria infections in the previous 3 months</td></tr><tr><td>pred_trim</td><td>Quarterly average number of mosquitoes child is exposed to</td></tr><tr><td>nutri_trim</td><td>Quarterly average nutrition scores</td></tr></tbody></table>

child, derived from a statistical predictive model based on climatic, entomological parameters, and characteristics of children's immediate surroundings as reported by [37].

Concerning the antibody quantification, two recombinant *P. falciparum* antigens were used to perform IgG subclass (IgG1 and IgG3) antibody quantification by Enzyme-Linked ImmunoSorbent Assay (ELISA) standard methods developed for evaluating malaria vaccines by the African Malaria Network Trust (AMANET [www.amanet148trust.org]). Protocol was described in detail [38].

**Data analysis** For our analysis, we use some of the data and we rename the proteins used in the study described above, for confidentiality reasons (some important findings are yet to be published). Thus, the proteins we use here, are named A1, A2, B and C, and are related to the antigens IgG1 and IgG3 as mentioned above in the description of the study. A1 and A2 are different domains of the same protein A, and C and D are two different proteins. Information contained in the multivariate longitudinal dataset of malaria are described in the Table 2.8, where Y denotes a protein which is one of the following:

$$ \text{IgG1\_A1, IgG3\_A1, IgG1\_A2, IgG3\_A2, IgG1\_B, IgG3\_B, IgG1\_C, IgG3\_C} \quad (2.43) $$

The aim of the analysis of these data is to evaluate the effect of the malaria infection on the child's immune (against malaria). Since the antigens which characterize the child's immune status interact together in the human body, we analyze the characteristics of the joint distribution of these antigens, conditional on the malaria infection and other factors of interest. The dependent variables are then provided by conc.Y (Table 2.8) which describes the level of the protein Y in the children at 3, 6, 9, 12, 15 and 18 months. All other variables in the Table 2.8 are covariates. We then have eight dependent variables which describe the longitudinal profile (in the child) of the proteins listed in Equation 3.48.

In the models that we fit to these data, we specify one random intercept by child and one random slope by child in the direction of the malaria infection. The illustration we do here is to jointly analyze each of the 28 pairs of proteins, in order to investigate if some profiles of proteins are independent, conditional on the configuration of the fitted model. After performing the bivariate
correlation test on all 28 bivariate models, the obtained p-values, with a Bonferroni correction, range from $4.16 \times 10^{-33}$ to 0.932. The p-value 0.932 is the only one which is not significant. This p-value corresponds to the pair of proteins (IgG3_A1, IgG1_B).

To investigate the general configuration of these proteins, in terms of correlations, we build their hierarchical cluster tree using $-\log(p-value)$ as dissimilarity. This hierarchical cluster tree is presented by the Fig 2.8.

The branch related to the IgG1 is different from the branch related to the IgG3. In other words, IgG1_A1, IgG1_A2, IgG1_B and IgG1_C are on the same branch which is different from the branch containing IgG3_A1, IgG3_A2, IgG3_B and IgG3_C (Fig 2.8). Relatively to both IgG1 and IgG3, A1 and A2 go together, and B and C also go together. These results are biologically very consistent, since A1 and A2 are domains of the same protein, and B and C are two different proteins. On the cluster (Fig 2.8), it also appears that the proteins IgG3_A1 and IgG1_B which are not significantly correlated (according to our bivariate test) are distant. Statistically, the model which may be used to jointly analyze these 8 protein profiles is not probably the model which contains all the 27 significant correlations, avoiding overfitting problems. Based on the results provided by the bivariate correlation test, it may be useful to perform a regularization procedure in the fitting of the full eight-variate model.

## 2.4 Conclusion

In the context of the multivariate linear mixed-effects model, we have suggested the more general expressions of the EM-based estimators than those used in the literature to analyze multivariate longitudinal data. These estimators fit the framework of the multivariate multilevel data analysis which, obviously, englobes the multivariate longitudinal data analysis framework. We also have built a likelihood ratio test based on these EM estimators to test the independence of two dimensions of the model. Furthermore, the simulation studies have validated the power of this test and have shown that this is an extremely sensitive test. In the context of longitudinal data, it allows to detect a modest correlation signal with a very small sample ($\rho = 0.3$, AUC= 0.81, with $n = 60$ subjects and $N = 600$ observations). In the simulation studies, the empirical distribution of the likelihood ratio statistic fits the $\chi^2(4)$. The asymptotic properties of likelihood ratio statistics, under nonstandard conditions, have been shown by [31] and [171]. These works have been generalized by [196] to cover a large class of estimation problems which allow sampling from non identically distributed random variables. The asymptotic distribution of the LR statistic derived by [196] is a mixture of chi-squared distributions. In the context of likelihood ratio tests for variance components in linear mixed-effects models, [74] used the results of [196] to prove that the proposed mixture of chi-squared distributions is the actual asymptotic distribution of such LR used as test statistics for null variance components with one or two random effects. Based on these works, Further theoretical investigations may be done to properly find out the asymptotic distribution of the likelihood ratio statistic in the case of this bivariate correlation test. Finally, we have illustrated the usefulness of the test on two different real-life data. The first dataset, which is of multivariate multilevel type, concerns the effects of school and classroom characteristics on pupils' progress in Dutch language and arithmetics, where the scores in language and arithmetics are the two response
variables which have been considered. Our method has yielded results that are consistent both with information in existing publications and with a conceptual understanding of the phenomenon. On this dataset, we have highlighted a joint effect between the scores in arithmetics and language within schools in the Netherlands. The second dataset, which is of longitudinal multivariate type, concerns a study of the effect of the malaria infection on the child's immune response in Benin. By jointly analyzing all the pairs of protein profiles of interest, we have plotted a hierarchical cluster tree of these proteins, using the bivariate correlation test. Information contained in this hierarchical cluster tree is consistent with the biological literature related to this issue.

The model as it is written is easily extendable to more dimensions despite a sparsity problem in choosing the parameterization of the covariance matrix or the precision matrix. Probably we could use this two-dimensional dependence test to structure a larger covariance matrix. The bivariate correlation test can help to construct iteratively, using a stepwise procedure, a parsimonious joint model containing all the components of $y$. This stepwise procedure may consist in adding to the constructing model, at each step, the significant correlation between two dependent variables. Using a model selection strategy, the model which fits more to the data will be retained. It could possibly be advantageous to turn to graphical LASSO type approaches to make a penalized estimation of this covariance (or precision) matrix. We could also resort to the rapid optimization methods such as that implemented in the lme4 [15] package, given the slow pace of the EM algorithm. It would be useful to assess the interest of this method compared to some heuristics such as the one which consists in setting one marginal response variable as a covariate of the other(s).

## Acknowledgements

We warmly thank the SCAC (Service de Coopération et d'Actions Culturelles) of the France Embassy in Benin, as well as the IRD (Institut de Recherche pour le Développement) for their financial support in the realization of this work.
<figure><img src="image_9.png" /><figcaption>Figure 2.7: Posterior distributions of random intercepts conditional on the data related to School 47 in the education dataset. Left panels assume independence across the two dimensions while right panels assume dependence. Top panels for the joint distribution of the random intercepts, middle panels for the joint distribution of random intercept in first dimension and random slope in the second dimension, bottom panels for the joint distribution of the random slopes.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 2.8: Hierarchical cluster tree on malaria-related proteins</figcaption></figure>
# PROFILED DEVIANCE FOR THE MULTIVARIATE LINEAR MIXED-EFFECTS MODEL FITTING

Dans ce chapitre, nous avons généralisé l'approche utilisée par Bates et al. (2014) dans le package lme4 du logiciel R pour estimer les paramètres d'un modèle linéaire multidimensionnel à effets mixtes. Cette généralisation est relative au cas où les résidus marginaux sont homoscédastiques avec la possibilité d'avoir des variances différentes. Les critères ML et REML sont explicitement écrits, mais pour des raisons de rapidité algorithmique, nous avons plutôt cherché à minimiser la déviance relativement aux paramètres de variance, avec une mise à jour itérative du vecteur des effets fixes jusqu'à convergence. Une étude de simulations a été menée et il en ressort que les estimations obtenues sont consistantes pour tous les paramètres du modèle. Ainsi la présente méthode surclasse l'estimation par l'algorithme EM, surtout en ce qui concerne les paramètres de variance. Ce travail a fait l'objet d'un article de recherche soumis dans une revue scientifique.
# Profiled deviance for the multivariate linear mixed-effects model fitting

Eric Houngla Adjakossa¹,²*, Mahouton Norbert Hounkonnou², Gregory Nuel¹

1 Laboratoire de Probabilités et Modèles Aléatoires /Université Pierre et Marie Curie, Case courrier 188 - 4, Place Jussieu 75252 Paris cedex 05 France

2 University of Abomey-Calavi, 072 B.P. 50 Cotonou, Republic of Benin

* ericadjakossah@gmail.com

## abstract

This paper focuses on the multivariate linear mixed-effects model, including all the correlations between the random effects when the marginal residual terms are assumed uncorrelated and homogeneous with possibly different standard deviations. The random effects covariance matrix is Cholesky factorized to directly estimate the variance components of these random effects. This strategy enables a consistent estimate of the random effects covariance matrix which, generally, has a poor estimate when it is grossly (or directly) estimated, using the estimating methods such as the EM algorithm. By using simulated data sets, we compare the estimates based on the present method with the EM algorithm-based estimates. We provide an illustration by using the real-life data concerning the study of the child's immune against malaria in Benin (West Africa).

## 3.1 Introduction

Linear mixed-effects model [62, 86, 90, 105, 192] has become a popular tool for analyzing univariate multilevel data which arise in many areas (biology, medicine, economy, etc), due to its flexibility to model the correlation contained in these data, and the availability of reliable and efficient software packages for fitting it [16, 82, 116, 141]. Univariate multilevel data are referred to as observations (or measurements) of a single variable of interest on several levels (school in a village which, in turn, is in a town), while multivariate multilevel data are characterized by multiple variables of interest measured at multiple levels. Examples include exam or test scores recorded for students across time, and multiple items at a single occasion for students in more than one school. Multivariate extension of the (single response variable-based) linear mixed-effects model is, indeed, having increasing popularity as flexible tool for the analysis of multivariate multilevel data [97, 159, 161, 197].

For the linear mixed-effects model, many methods for obtaining the estimates of the fixed and the random effects have been proposed in the literature. These methods include Henderson's mixed model equations [91], approaches proposed by [75] as well as techniques based on two-stage regression, Bayes estimation, etc. For details, see [168, Section 7.4c] and [155]. Concerning
the variance parameters estimation in linear mixed-effects model, the discussed methods in the literature include the ANOVA method for balanced data which uses the expected mean squares approach [169, 170]. For unbalanced data, [149] proposed the minimum norm quadratic estimation (MINQUE) method, where the resulting estimates are translation invariant under unbiased quadratic forms of the observations. [109] gave another method of estimating variance parameters using extended quasi-likelihood, i.e. gamma-log generalized linear models. For more details on these parameters' estimation methods in the linear mixed-effects model, see the paper of [81]. Beside all the methods cited earlier, come the Maximum Likelihood (ML) and the Restricted Maximum Likelihood (REML) methods. ML and REML methods are the most popular estimation methods in the linear mixed-effects model [114]. The main attraction of these methods is that they can handle a much wider class of variance models than simple variance components [81].

In the multivariate linear mixed-effects model, ML and REML estimates are frequently approached through iterative schemes such as EM algorithm [6, 44, 125, 161, 172]. This avoid the difficulties related to the direct calculating of the parameters' likelihood, since the random effects are not observed, without ignoring the flexible computationally of these algorithms. Despite the existence of valid theorems which show the asymptotic convergence of the sequences produced by these algorithms toward ML estimates [44], in practice this may not always work exactly as expected.

In this paper, we focus on the multivariate linear mixed-effects model, including all the correlations between the random effects while the marginal residuals are assumed independent homogeneous with possibly different standard deviation. The class of multivariate mixed-effects models considered here assumes that the random effects and the residuals follow Gaussian distributions, and the dependent variables are continuous. In this model, our approach consists in directly calculating the likelihood of the model's parameters. This likelihood is used to obtain the ML estimates or the REML estimates through the provided REML criterion. This strategy may explain the high quality of the estimates of both fixed effects parameters and random effects' variance parameters as well as residual variance parameters. This approach may be viewed as a generalization of the approach proposed by [16] under the R software [147] package named lme4.

## 3.2 Multivariate linear mixed-effects model

For the sake of simplicity we focus on the bivariate case ($d = 2$) in most of the paper, but the generalization to higher dimensions ($d > 2$) is straightforward. Thus, in dimension 2, the model is the following:

$$ 
\begin{aligned}
y_1 &= X_1\beta_1 + Z_1\gamma_1 + \varepsilon_1, \\
y_2 &= X_2\beta_2 + Z_2\gamma_2 + \varepsilon_2,
\end{aligned} 
\qquad (3.1) 
$$
where

$$ \gamma = \begin{pmatrix} \gamma_1 \\ \gamma_2 \end{pmatrix} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{\Gamma} = \begin{pmatrix} \Gamma_1 & \Gamma_{12} \\ \Gamma_{12}^\top & \Gamma_2 \end{pmatrix} \right), \quad \varepsilon = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \end{pmatrix} \sim \left( \mathbf{0}, \begin{pmatrix} \sigma_1^2 \mathbf{I}_N & 0 \\ 0 & \sigma_2^2 \mathbf{I}_N \end{pmatrix} \right). \quad (3.2) $$

For the sake of simplicity, we write $\gamma \sim \mathcal{N}(\mathbf{0}, \mathbf{\Gamma})$ to mean that $\gamma$ is a realization of a random vector which is $\mathcal{N}(\mathbf{0}, \mathbf{\Gamma})$ distributed. For $k \in \{1, 2\}$, $\beta_k$ and $\gamma_k$ denote respectively the fixed effects and the random effects vector of covariates, while $\varepsilon_k$ is the marginal residual component in the dimension $k$ of the model. $X_k$ is a matrix of covariates and $Z_k$ a covariates-based matrix of design. $\dim(X_k) = N \times p_k$ and $\dim(Z_k) = N \times q_k$, where $N$ is the total number of observations. $p_k$ and $q_k$ are, respectively, the number of fixed effect related covariates and the number of random effect related covariates in the dimension $k$ of the model. $\mathbf{y} = (y_1^\top, y_2^\top)^\top$ is the vector of marginal observed response variables of the model. We assume that $\mathbf{y}$ is a realization of a random vector $\mathcal{Y}$ and belongs to $\mathbb{R}^{2N}$. The bold symbols represent parameters, or vectors, of multiple dimensions (i.e. $\Gamma_1$ concerns dimension 1 of the model while $\Gamma$ concerns both dimensions).

$\Gamma_1$ and $\Gamma_2$ are the variance-covariance matrices of $\gamma_1$ and $\gamma_2$, respectively. $\Gamma_1$ and $\Gamma_2$ must be, indeed, positive semidefinite. It is then convenient to express the model in terms of the relative covariance factors, $\Lambda_{\theta_1}$ and $\Lambda_{\theta_2}$, which are $q_1 \times q_1$ and $q_2 \times q_2$ matrices, respectively. $\Lambda_{\theta_1}$ is a block diagonal matrix. Each element in the diagonal of $\Lambda_{\theta_1}$ is a lower triangular matrix whose nonzero entries are the components of the vector $\theta_1$. That is, $\theta_1$ generates the symmetric $q_1 \times q_1$ variance-covariance matrix $\Gamma_1$, according to

$$ \Gamma_1 = \sigma_1^2 \Lambda_{\theta_1} \Lambda_{\theta_1}^\top. \quad (3.3) $$

Same as $\theta_2$ which generates $\Gamma_2$ according to

$$ \Gamma_2 = \sigma_2^2 \Lambda_{\theta_2} \Lambda_{\theta_2}^\top. \quad (3.4) $$

In Equations 3.3 and 3.4, $\sigma_1^2$ and $\sigma_2^2$ are the same marginal residual variances used in the model expression (see Equation 3.2). Using the variance-component parameters, $\theta_1$ and $\theta_2$, the marginal random effects, $\gamma_1$ and $\gamma_2$, are expressed as

$$ \gamma_1 = \Lambda_{\theta_1} u_1, \quad \gamma_2 = \Lambda_{\theta_2} u_2, \quad (3.5) $$

such that

$$ \mathbf{u} = \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} \sim \mathcal{N}(\mathbf{0}, \Sigma_{\mathbf{u}}), \quad \text{with} \quad \Sigma_{\mathbf{u}} = \begin{pmatrix} \sigma_1^2 \mathbf{I}_{q_1} & \sigma_1 \sigma_2 \rho \\ \sigma_1 \sigma_2 \rho^\top & \sigma_2^2 \mathbf{I}_{q_2} \end{pmatrix}. \quad (3.6) $$