Formation SFP 2018 : les modèles linéraires mixtes

1. Prérequis

1.1. La variance et la covariance

Les notions de variance et de covariance sont considérées être connues. Pour rappel, voici les formules de la variance et de la covariance et la manière d’obtenir ces informations en utilisant R.

Variance	Covariance
\(s^2= \sum\frac{(x_i-\mu)^2}{N-1}\)	\(cov_xy= \sum\frac{(x_i-\mu_x) (y_i-\mu_y)}{N-1}\)

var(bfi[,1:4], na.rm=T)->ex.var
ex.var

##            A1         A2         A3         A4
## A1  1.9747487 -0.5651005 -0.4861049 -0.3050179
## A2 -0.5651005  1.3800668  0.7434176  0.5822290
## A3 -0.4861049  0.7434176  1.7001245  0.7013222
## A4 -0.3050179  0.5822290  0.7013222  2.2054776

Dans cette matrice, la diagonale représente les variances des variables.

diag(ex.var) # les variances

##       A1       A2       A3       A4 
## 1.974749 1.380067 1.700125 2.205478

Ce qui est en dehors de la diagonale représente les covariances :

ex.var[lower.tri(ex.var, diag=T)]<-NA
ex.var # les covariances

##    A1         A2         A3         A4
## A1 NA -0.5651005 -0.4861049 -0.3050179
## A2 NA         NA  0.7434176  0.5822290
## A3 NA         NA         NA  0.7013222
## A4 NA         NA         NA         NA

1.2. La distribution normale

Une distribution suit une distribution normale quand la distribution d’une variable numérique suit une distribution en cloche comme ci-desssous.

Elle se caractérise par la notation suivante :

\[ X \sim N(\mu,\sigma^2) \] La notation signifie que la variable aléatoire X suit une distribution normale qui dépend de deux paramètres : la moyenne et la variance.

2. Le modèle linéaire mixte (MLM)

2.1. Quelques généralités

2.1.1. Bref historique des logiciels permettant de réaliser des modèles linéaires mixtes.

• 1992 : SAS est le premier logiciel où il est possible de faire des modèles linéaires mixtes.
• 1998 : Le package ‘nlme’ de R apparaît.
• 2001 : SPSS permet de faire des modèles linéaires mixtes
• 2005 : STATA sort une version avec un module où il est possible de faire les modèles linéaires mixtes et le package ‘lme4’ est mis à disposition

2.1.2. Définition

Un modèle linéaire mixte est un modèle pour lequel le modèle comprend à la fois des effets fixes et des effets aléatoires. Les MLM incluent des variables à effets fixes et aléatoires. Le mélange entre les deux types de facteurs dans un même modèle est à l’origine du nom. Les effets fixes décrivent les relations entre les covariables et la variables dépendante pour une population entière, les effets aléatoires sont spécifiques à l’échantillon.

En d’autres termes, - un effet aléatoire est effet dont nous ne voulons pas généraliser les propriétés (les modalités ont été choisies de manière aléatoire dans quelque chose de plus grand).

• un effet fixe est un effet dont on veut généraliser les propriétés. Il s’agit de la variable manipulée. Les niveaux de cette variable ont été choisi de manière spécifique.

Il est important de comprendre qu’une variable peut être considérée comme un effet fixe ou un effet aléatoire en fonction de l’hypothèse qui va être testée.

IMPORTANT 1 : les effets aléatoires doivent nécessairement être des variables catégorielles. L’utilisation de nombres peut créer de l’ambiguïté. Il est donc préférable de s’assurer que le format de la variable aléatoire soit une variable catégorielle.

IMPORTANT 2 : six modalités de cette variable catégorielle semble être un minimum pour pouvoir estimer de manière pertinente une variance.

2.1.2.1 les facteurs fixes.

Les facteurs fixes sont les variable pour lesquelles tous les niveaux (i.e., toutes les conditions) qui sont d’intérêt sont inclus dans l’étude. Les variables peuvent être qualitatives comme le genre, ou des variables d’échantillonnage de l’étude comme la région. Les niveaux sont choisis de sorte à représenter des conditions spécifiques qui peuvent être utilisées pour définir des contrastes.

2.1.2.2 les facteurs aléatoires

Les effets aléatoires sont échantillonnés de manière aléatoire dans l’ensemble des niveaux possibles de la population étudiée. On n’a pas tous les niveaux de la population, mais l’objectif est de pouvoir inférer les propriétés de l’échantillon à une population entière.

Contrairement aux facteurs fixes, les niveaux de la variables aléatoire ne représentent pas des conditions choisies de manière spécifique pour répondre aux objectifs de l’étude.

2.1.2.3 les effets fixes vs les effets aléatoires

Les effets fixes sont représentés par des coefficients de régression. Ces effets décrivent les relations entre la variable dépendante et les prédicteurs (des facteurs fixes ou des covariables continues). Ces effets permet d’identifier, par exemple, les différences de moyennes entre des groupes sur la VD pour les facteurs fixes, ou le lien entre une covariable continue et la variable dépendante. On fait l’hyptohèse que les effets fixes sont inconnus et que nous les estimons sur la base des données récoltées.

Les effets aléatoire sont spécifique à un niveau donné du facteur aléatoire. Ces effets représentent une déviation de la relation décrite par les effets fixes.

2.1.3. Modèle linéaire mixte, modèle multiniveau, analyse de courbe de croissance.

On parle d’analyse de courbe de croissance lorsque le modèle linéaire mixte a pour objectif d’évaluer la trajectoire de la variable dépendante différenciée à la fois en fonction des personnes et en fonction du temps.

Les modèles mixtes sont aussi appelés modèles multiniveaux lorsqu’une des variables est organisée en cluster : un groupe de participant est choisi à l’intérieur d’une variable de plus haut niveau. Un exemple classique d’un modèle multiniveau est celui des enfants qui appartiennent à une classe. Ils vont plus se ressembler car ils ont eu le même enseignant, les enfants issus d’une même école vont plus se ressembler en raison par exemple de la culture d’établissement, de la localisation de l’école…

3. Mise en application avec R.

3.1. Les packages R pour les MLM.

Nous allons utiliser essentiellement deux packages R qui permettent d’estimer les modèles linéaires mixtes : le package ‘nlme’ et ‘lme4’

Ces deux packages ont été écrits par Douglas M. Bates.

Chacun de ces deux packages ont des avantages et des inconvénients.

Package nlme	Package lme4
Plus lent	Plus rapide
Permet de tester l’hétéroscédasticité	Ne permet pas de tester les variances
Permet de tester les covariances	Ne permet pas de modéliser les covariances
Plus difficile de modéliser plusieurs facteurs aléatoires	Permet de modéliser plusieurs facteurs aléatoires
Ne permet pas de tester les modèles non linéaires	Possibilité de modéliser des facteurs non linéaires
Ne corrige pas les ddl dans la table des F et t	Corrige les ddl pour les F et les t

3.2. Débuter avec R : problèmes et solutions

R un est logiciel de programmation et il est souvent difficile pour les débutants de s’y mettre car ils n’ont pas de notion de programmation. Ainsi, des erreurs faciles à éviter produisent des erreurs qu’il est souvent difficile à comprendre quand on commence. Pour les éviter, il est importer d’adopter des règles de bonnes pratiques.

3.2.1. Garder une trace : utiliser un script

Puisque nous allons utiliser R en ligne de commande, on peut utiliser R ou Rstudio. La différence fondamentale entre R et Rstudio est le fait que Rstudio est plus convivial dans la mesure où un code couleur va apparaître. Par ailleurs, vous pouvez importer assez aisément des données dans Rstudio grâce à l’onglet “import Dataset”. Il est donc préférable d’utiliser un Rstudio plutôt que R (bien que les deux soient en réalité strictement identiques).

Une fois Rstudio (ou R) ouvert, il est préférable d’utiliser un script et d’éviter de taper directement les lignes de commandes dans la consoles pour deux raisons : - il est plus difficile de corriger une erreur dans la console - le script vous permettra de garder une trace de ce que vous faites.

En effet, il est possible de commenter à l’aide de dièse (#) les lignes de commandes que vous allez utiliser. Par exemple :

# Le fait d'utiliser un dièse me permet de commenter
getwd() # permet d'obtenir le répertoire de travail

## [1] "C:/Users/Nicolas/Documents/reims/recherches/SFP2018/LME"

En résumé, Utiliser un script dans R plutôt que la console permet de corriger plus facilement les erreurs, de sauvegarder les lignes de commandes utilisées, et d’introduire des commentaires (en précédant la ligne par le sigle #).

Exercice 1 : lancez Rstudio, et ouvrez le script LME.SFP2018.R

3.2.2. Charger les packages dont on a besoin

R est un logiciel avec un coeur qui a relativement peu de fonctionnalités. En revanche, il existe un très vaste ensemble de packs de fonctions supplémentaires qui ont été développés par la communauté. Ces packs s’appellent des packages. Une des raisons les plus fréquentes des bugs est qu’un (ou plusieurs packages) n’ont pas été chargé. Pour éviter ce problème, le mieux est de charger l’ensemble des packages dont on va avoir besoin dès le début de la session.

Il existe deux moyens pour charger un package : - la fonction “library” - la fonction “require” Par exemple, on peut charger le package ‘nlme’ de la manière suivante :

require(nlme)

## Loading required package: nlme

Ou de la manière suivante :

library(nlme)

La différence pratique entre les deux est que “require” ne renvoie pas de message d’erreur quand le package n’a pas été chargé correctement contrairement à “library”. Ainsi, il est préférable d’utiliser “library”. Cependant, les packages ne sont pas installés d’emblée sur les ordinateurs. Il est donc nécessaire de les installer au préalable. Pour cela, on utilise la fonction install.packages avec le nom du package entre guillemets. Par exemple, on installe le package “nlme” de la manière suivante :

install.packages("nlme", repos = "http://mirror.ibcp.fr/pub/CRAN" ) # l'argument repos est accessoire 

## Installing package into 'C:/Users/Nicolas/Documents/R/win-library/3.6'
## (as 'lib' is unspecified)

## Warning: package 'nlme' is in use and will not be installed

Exercice 2 : dans le script qui vous est fourni, identifiez et lancez les lignes de commandes qui permettent d’installer et de lancer l’ensemble des packages dont nous aurons besoin dans le cadre de cette formation.

3.2.3. L’importation des données

Quand vous voulez importer un fichier de données, il arrive que R ne trouve pas le fichier ou vous indique que le chemin d’accès est inexistant/introuvable … Le problème est que le chemin d’accès a été mal spécifié.

Pour éviter ce problème : - il faut commencer par spécifier le répertoire de travail en choisissant le dossier dans lequel se trouve votre base de données. Pour connaître le répertoire de travail, on utilise la fonction suivante :

getwd()

## [1] "C:/Users/Nicolas/Documents/reims/recherches/SFP2018/LME"

et celle-ci pour spécifier le répertoire de travail :

setwd("C:/Users/Nicolas/Documents/reims/recherches/SFP2018/LME")

Lorsque vous fixer le répertoire de travail, veillez à respecter la casse (MAJUSCULE/minuscule).

Une fois dans le répertoire de travail, il faut s’intéresser au nom de la base de données. Plus spécifiquement, il faut bannir : * les accents * Les espaces * Les caractères spéciaux * Les noms trop longs

Ces règles valent également pour le nom des variables dans la base de données. Le nom d’une variable n’est pas le chapitre d’un ouvrage. Il ne doit pas faire plus d’une trentaine de caractères.

Exercice 3: en utilisant l’onglet ‘import Dataset’ (voir Figure 1). de Rstudio, importez les données contenues dans la feuille “rat” du fichier excel “SFP2018.xls” et donnez-lui le nom de “rat”.

Figure 1. Onglet Rstudio permettant d’importer un jeu de données.

3.2.3. Stocker des informations dans la mémoire de R

Généralement, les sorties d’une fonctions ne sont pas directement utilisables. Il est encore nécessaire d’exécuter d’autres fonctions sur les résultats renvoyés par une fonction. Il faut donc stocker l’information dans la mémoire de R. Cela se fait aisément par le symbole “<-” ou “->”. Par exemple, on va stocker le chiffre 5 dans un objet appelé “exemple” de la manière suivante :

exemple<- 5 # les espaces ne sont pas importants. Est strictement identique à 5->exemple
exemple # faire attention à l'orthographe et la casse. "Exemple" renvoie un message d'erreur

## [1] 5

Il est facile de savoir ce qu’il y a dans la mémoire de R grâce à Rstudio. En effet, ce qui est stocké dans la mémoire de R apparaît en haut à droite de la fenêtre. Comme dans l’exemple présenté dans la figure 2 où seule un jeu de données appelé “rat” est présent dans la mémoire de R.

Figure 2. Liste des objects stockés dans la mémoire de R. Seul un jeu de données appelé “rat” est présent.

Ainsi, un des points souvent oubliés est de stocker la sortie de résultats de R dans la mémoire pour appliquer d’autres fonctions par la suite à cette sortie de résultats.

4. La formalisation du modèle

4.1. Description de l’exemple de travail numéro 1 : les portées de rats

Les données de la feuille “rat” du fichier excel “sfp2018.xlsx” proviennent d’une étude dans laquelle 30 rattes femelles ont été assignés de manière aléatoire à une des trois doses de traitement (haut, bas, contrôle). Le but de l’étude est de comparer le poids à la naissance des rats issus de ces trois groupes. Bien que, au départ, il y avait 10 rats dans chacun des groupes, 2 rats sont décédés durant l’étude. De plus, les portées variaient considérablement, allant de 2 à 18 rats. Pour ces raisons, le plan de l’étude n’était pas équilibré.

library(readxl)
rat <- read_excel("formation.LME.xlsx", sheet="rat")
str(rat)

## Classes 'tbl_df', 'tbl' and 'data.frame':    322 obs. of  6 variables:
##  $ id       : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ weight   : num  6.6 7.4 7.15 7.24 7.1 6.04 6.98 7.05 6.95 6.29 ...
##  $ sex      : chr  "Male" "Male" "Male" "Male" ...
##  $ Litter   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Lsize    : num  12 12 12 12 12 12 12 12 12 12 ...
##  $ Treatment: chr  "Control" "Control" "Control" "Control" ...

Dans un jeu de données, les variables “integer” sont des nombres entiers, des variables “num” sont des variables numériques, des variables “chr” sont des chaînes textuelles et des variables “factor” sont des variables catégorielles. Les variables dans le jeu de données sont décrites dans le Tableau 1.

Tableau 1. Description du jeu de données ‘rat’.

Variable	Description
id	Identifiant du rat. Il y a 322 rats en tout. Doit être catégorielle
weight	Poids du rat. Variable numérique
sex	Sexe du rat. Variable catégorielle
Litter	Portée à laquelle appartient le rat. Doit être catégorielle
Lsize	Taille de la portée. Variable numérique
Treatment	Traitement auquel la mère a été soumis. Variable catégorielle

Il apparaît que la structure des données ne correspondent pas tout à fait, en particulier pour la variable “id” et la variable “Litter” qui sont nnumériques au lieu d’être catégorielle. On va donc modifier la nature de ces variables. En utilisant la fonction “factor”. Il faut donc dire à R que la variable “id” dans la base de données “rat” est un facteur et le stocker dans la mémoire de R. Voici comment procéder :

rat$id<-factor(rat$id) # le symbole dollar permet d'indiquer la variable "id" est dans le jeu de données "rat"

Exercice 4: transformez les variables qui doivent l’être en variables catégorielles

On peut voir les premières lignes de la base de données grâce à la fonction “head”.

head(rat)

## # A tibble: 6 x 6
##   id    weight sex   Litter Lsize Treatment
##   <fct>  <dbl> <fct> <fct>  <dbl> <fct>    
## 1 1       6.6  Male  1         12 Control  
## 2 2       7.4  Male  1         12 Control  
## 3 3       7.15 Male  1         12 Control  
## 4 4       7.24 Male  1         12 Control  
## 5 5       7.1  Male  1         12 Control  
## 6 6       6.04 Male  1         12 Control

4.2. Comparaison entre le modèle linéaire classique et le modèle linéaire mixte

4.2.1 L’anova vs le MLM : similitudes

Les modèles linéaires mixtes sont des alternatives au modèle linéaire, et en particulier aux anovas. Tout comme les anovas, les modèles linéaires mixtes permettent de déterminer si des variables catégorielles ont un impact sur une variable dépendante continue.

Tout comme les anovas, le résidu doit être distribué en suivant un distribution normale.

Les prédicteurs peuvent être de nature catégorielle, ordinale ou continue.

Le modèle linéaire classique est défini de la manière suivante :

\[ Y=β_0 +β_1 X+ ε \]

Où Y représente les mesures observées sur la variable dépendante, \(\beta_0\) est l’intercept et \(\beta_1\) est le coefficient de régression. Dans ce cas, \(\epsilon\) représente le résidu, c’est-à-dire ce que nous ne sommes pas en mesure de modéliser.

Si nous nous intéressons à un individu i, la prédiction pour cet individu i est

\[ Y_i=β_0 +β_1 X_i+ ε \] Cela signifie que, pour un individu i, la valeur qu’on peut prédire pour la variable dépendante correspond à l’intercept auquel on ajoute la valeur de l’individu i pour la variable indépendante multiplié par le coefficient de régression. Entre cette valeur prédite de manière théorique et la valeur réellement obserbée, il y a une différence. La différence entre la valeur prédite et la valeur observée est le résidu et ce résidu n’est pas modélisé de manière explicite dans le modèle linéaire.

Dans R, on modélise le modèle linéaire par la fonction “lm”, l’acronyme de linear model. Un modèle est défini par selon la structure suivante : variable dépendante ~ prédicteur 1 + prédicteur 2.

Donc, on peut définir un modèle linéaire où le poids des rats est défini par le traitement de la manière suivante :

lm.out<- lm(weight~ Treatment,data=rat)

Il faut remarquer que la sortie de la fonction “lm” a été stocké dans un objet qu’on va appelé “lm.out”. Cela permet d’obtenir des résultats complémentaires. En particulier, on peut obtenir de la table de l’anova avec la fonction “anova” :

 anova(lm.out)

## Analysis of Variance Table
## 
## Response: weight
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## Treatment   2  13.202  6.6011  17.353 7.018e-08 ***
## Residuals 319 121.349  0.3804                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On peut également obtenir les coefficients de la droite de régression avec la fonction “summary”.

 summary(lm.out)

## 
## Call:
## lm(formula = weight ~ Treatment, data = rat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.64473 -0.34743 -0.03694  0.31007  2.00527 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.32473    0.05389 117.370  < 2e-16 ***
## TreatmentHigh -0.43919    0.09357  -4.694 3.99e-06 ***
## TreatmentLow  -0.39640    0.07696  -5.151 4.55e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6168 on 319 degrees of freedom
## Multiple R-squared:  0.09812,    Adjusted R-squared:  0.09247 
## F-statistic: 17.35 on 2 and 319 DF,  p-value: 7.018e-08

Ainsi, notre droite de régression devient : \[ Poids_i= 6.32+ (-0.44) * Traitement_{.High} + (-0.40) * Traitement_{.low} \]

Plusieurs points méritent un peu d’attention :

1. R considèrent par défaut les modalités d’une variable catégorielle par ordre alphabétique. Ainsi, la ligne de base (la condition de référence) est le groupe de contrôle dans notre exemple. En d’autres termes, l’équation ci-dessus signifie que le groupe de contrôle pèse en moyenne 6.32 grammes, que les rats dont la mère a reçu une dose élevée pèse en moyenne 0.44 grammes en moins que le groupe de contrôle, et un rat dont la mère a reçu une dose faible pèse en moyenne 0.40 grammes en moins que le groupe de contrôle.

2. Le nombre de coefficients de régression (en plus de l’intercept) correspond au nombre de modalités -1

3. Les coefficients de régressions correspondent à une modalité spécifique d’une variable catégorielle.

4. ce type de contrastes correspond à ce qui est appelé dans R à “traitement” et correspond à la structure suivante.

contrasts(rat$Treatment)

##         High Low
## Control    0   0
## High       1   0
## Low        0   1

5. Par défaut, le niveau de référence dans R est la première modalité alphabétique (ou numérique).

Dans un modèle linéaire mixte, contrairement au modèle linéaire classique, on va formaliser de manière explicite la variable aléatoire.

\[ Y_i=β_0+β_1 X_i+ u_i Z+ϵ_i \]

Ainsi, par rapport au modèle linéaire, l’équation du modèle a été enrichie de \(u_i Z\) où \(u_i\) représente le “coefficient de régression de la variable aléatoire” pour l’individu i et Z étant la valeur de l’individu i.

De manière concrète, on va modéliser le MLM avec la fonction ‘lme’ du package ‘nlme’ ou ‘lmer’ du package ‘lme4’. Pour notre comparaison, nous utiliserons uniquement la fonction lme. On formalise un MLM de manière assez similaire au modèle linéaire. Cependant, il y a un argument supplémentaire, l’argument ‘random’ qui permet de préciser le facteur aléatoire. On le définit de la manière suivante : \(random=1|var.aléatoire\). Le chiffre 1 indique la variable aléatoire est liée à l’intercept de notre modèle (i.e., il y a un intercept pour chaque modalité de la variable aléatoire), la barre verticale indique qu’il s’agit d’un effet aléatoire et on donne le nom de la variable. Dans notre exemple, la variable aléatoire est l’identifiant du rat.

lme.out<-lme(weight~Treatment, data=rat, random=~1|id)  # identifier les effets fixes et les effets aléatoires. 

Que ce soit pour la table de l’anova,

anova(lme.out)

##             numDF denDF   F-value p-value
## (Intercept)     1   319 31300.817  <.0001
## Treatment       2   319    17.353  <.0001

ou pour les coefficients de régression,

summary(lme.out)

## Linear mixed-effects model fit by REML
##  Data: rat 
##        AIC      BIC    logLik
##   620.8482 639.6741 -305.4241
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept) Residual
## StdDev:   0.5774987 0.216562
## 
## Fixed effects: weight ~ Treatment 
##                   Value  Std.Error  DF   t-value p-value
## (Intercept)    6.324733 0.05388735 319 117.36954       0
## TreatmentHigh -0.439194 0.09357464 319  -4.69352       0
## TreatmentLow  -0.396399 0.07696054 319  -5.15069       0
##  Correlation: 
##               (Intr) TrtmnH
## TreatmentHigh -0.576       
## TreatmentLow  -0.700  0.403
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.50563322 -0.19779199 -0.02102742  0.17651948  1.14158862 
## 
## Number of Observations: 322
## Number of Groups: 322

on constate que les coefficients pour les effets fixes correspondent aux coefficients pour le modèle linéaire classique. En plus des informations fournies par le modèle classique, nous avons des indications sur la variance des effets aléatoires. Donc, juusqu’à présent, nous n’avons pas présenté pourquoi passer au modèle plus complexe du modèle linéaire mixte vu que ce sont les mêmes résultats. L’objectif était de montrer que, dans la pire des situations, le MLM est une anova classique.

Pour notre exemple de travail, le modèle que nous allons vouloir tester est un peu plus complexe : \[ Poids_ij = β0 + β1 × TRAIT1j + β2 × TRAIT2j + β3 × SEX1ij+ \]

\[ β4 × Taille.Portée + β5 × TRAIT1j × SEX1ij+ β6 × TRAIT2j × SEX1ij+uj + εij\]

où j représente la portée et i le bébé rat au sein de la portée.

Dans cet exemple, la portée va être considérée comme une variable aléatoire. On peut donc associer un intercept aléatoire, représentant la déviation du poids pour chaque portée par rapport à l’intercept global.

4.2.2. Le MLM vs l’anova : les différences

4.2.2.1. Les conditions d’application

En plus du résidu, les variables aléatoires doivent aussi suivre une distribution multinormale.

En revanche :

• les variances peuvent être hétérogènes (et on peut le modéliser).

• Les résidus peuvent ne pas être indépendants (il peut y avoir des corrélations sur les résidus).

• Il est possible d’ajuster le modèle avec une large sélection de structure de covariance parcimonieuse (ce qui est plus efficace que d’estimer la structure complète de variance-covariance).

• On peut modéliser la distribution de la variable dépendante (qui n’est pas nécessairement continue).

• La sphéricité de la matrice de covariance n’est pas requise pour les modèles ayant des facteurs en mesure répétée.

1. on peut suivre des individus dans le temps même si le nombre de mesures est différent (on n’enlève pas tout l’individu parce qu’on a une valeur manquante). Dans le modèle linéaire mixte, toutes les observations disponibles pour un individu sont utilisées dans l’analyse.

2. lorsqu’on utilise une anova, les niveaux pour le facteur temps doivent être les mêmes. Dans le modèle linéaire mixte, les mesures de temps peuvent varier entre les participants. En revanche, il faut tout de même que la plupart des informations soit disponible à chaque mesure de temps.

• Le plan ne doit pas être équilibré, même pour les facteurs à mesure répétée. En d’autres termes, il peut y avoir des valeurs manquantes sur les facteurs à mesure répétée sans perdre des participants.

• On peut modéliser des plans expérimentaux pour lesquels certains facteurs sont partiellement emboîtés.

En résumé, toutes ces différences entre le MLM et l’ANOVA permet de comprendre que les MLM sont bien plus flexibles que les anovas.

En revanche, les propriétés du modèle linéaire mixte concernant les valeurs manquantes ne sont possibles qu’à la condition que les valeurs soient manquante de manière aléatoire (MAR - missing at random).

Le MAR signifie que la probabilité d’avoir des valeurs manquantes sur une variable donnée puisse dépendre d’autres informations mais ne dépend pas de données qui auraient pu être observées mais sont en fait manquantes. Par exemple, si une personne ne reporte pas son poids parce qu’elle se trouve trop grosse, ce n’est pas un MAR. En revanche, si le jour où on enregistre le poids, la personne est malade, il s’agit d’un MAR.

Lorsqu’on utilise des données en cluster (les élèves qui sont emboîtés dans une classe), les clusters peuvent avoir des tailles inégales ou il peut y avoir des MAR. Cela entraîne des designs qui ne sont pas équilibrés. Cela ne pose pas de souci pour les MLM (tant que les données manquantes sont des MAR). On peut donc faire des inférences qui ont du sens même sur des jeu de données qui ne sont pas équilibrés.

Pour illustrer les propriétés des modèles linéaires mixtes, dans notre exemple, on comprend aisément que les unités d’analyses (les rats qui naissent au sein d’une même portée) vont présenter des résultats qui seront corrélés puisque les fœtus partageaient le même environnement maternel. Pourtant, dans chaque portée, il n’y a pas le même nombre de rats, et la répartition en fonction du sexe n’est pas la même d’une portée à l’autre.

Ainsi, dans le modèle, la portée va donc être un effet aléatoire (impliquant que les observations d’une même portée sont corrélées) et nous pourrons spécifier si nous le souhaitons les corrélations au sein des portées.

4.2.2.2. Les covariables

• Les modèles mixtes permettent d’avoir à la fois des covariables au niveau individuel (tel que l’âge ou le sexe) et au niveau du cluster (telle que la taille du cluster) lorsqu’on ajuste l’effet aléatoire.

• On peut avoir des covariables qui varient en fonction du temps dans le modèle.

4.2.3. Nomenclature dans les modèles linéaires mixtes : les variables croisées (crossed) et emboîtées (nested)

Une variable croisée est une variable pour laquelle les mesures sur la variable dépendante sont faites sur les mêmes individus pour les différentes modalités de la variable (mesures répétées).

La littérature sur les modèles linéaires mixtes distingue les plans en mesure répétées et les études longitudinales. La différence fondamentale entre les deux est que, dans une étude longitudinale, on mesure les données à plusieurs moments, généralement sur une période relativement longue. La conséquence est que la perte expérimentale est plus importante.

Si la distinction entre entre mesure répétée et étude longitudinale n’est pas critique, le point important à noter est que, dans les deux cas, la variable dépendante est mesurée plus qu’une fois pour chaque unité de mesure, et ces mesures sont vraisemblablement corrélées.

Une variable emboîtée est un variable pour laquelle les mesures sur la variable dépendante sont réalisées sur des individus différents entre les différentes modalités de la variable. On parle aussi dans ce cas de données en cluster.

4.3. Les fonctions R permettant d’ajuster un modèle linéaire mixte

Comme évoqué précédemment, deux fonctions permettent d’ajuster une MLM : la fonction ‘lme’ du package ‘nlme’ et la fonction ‘lmer’ du package ‘lme4’.

Nous allons décrire les deux fonctions et modéliser notre modèle lié au jeu de données ‘rat’ avec ces deux fonctions.

4.3.1. La fonction ‘lme’

Pour la fonction ‘lme’ du package ‘nlme’, nous nous concentrerons sur 7 arguments qui sont décrits dans le Tableau 2.

Tableau 2. Description des arguments de la fonction ‘lme’.

Argument	Description	Illustration
fixed	Modèle qu’on veut tester	lme(weight ~ Treatment + sex + Lsize +Treatment:sex,
data	Nom du jeu de données	data = rat,
random	Nom de l’effet aléatoire	random = ~ 1 |Litter ,
method	Choix de l’algorithme (“ML” ou “REML”)	method = “REML”,
na.action	Manière de gérer les valeurs manquantes	na.action=na.omit,
correlation S	pécification la structure de corrélation	correlation = NULL,
weights	Modélisation de la structure de variance	weights = varIdent(form = ~1 | Treatment))

Ce qu’il faut identifier sur la description des arguments de la fonction : - Dans R, les arguments d’une fonction sont séparés par une virgule. - Il n’est pas nécessaire de spécifier un argument qu’on n’utilise pas. Par exemple, en l’occurrence, il n’y a ni valeurs manquantes, ni spécification de la structure de corrélations. On pourrait donc écrire le modèle comme ceci :

Modele.lme<-lme(  weight ~ Treatment + sex + Lsize +Treatment:sex, data = rat, 
            random = ~ 1 |Litter ,   method = "REML", weights = varIdent(form = ~1 | Treatment)) 

• La méthode est l’algorithme de calcul. Il s’agit du maximum de vraisemblance (ML) ou du maximum de vraisemblance restreint (REML). On utilise le maximum de vraisemblance quand l’hypothèse porte sur un effet fixe et le maximum de vraisemblance restreint quand l’effet porte sur un effet aléatoire. Pour en savoir plus, veuillez vous reporter à la section consacrée au maximum de vraisemblance en supplément de cette formation.

• Dans notre exemple, nous avons voulu modéliser des variances hétérogènes sur la variabe traitement.

• Le facteur aléatoire de la portée (Litter) est explicitement identifié comme étant un facteur aléatoire (random)

• Quand on utilise un logiciel avec boîte de dialogue, l’usage habituel est de modéliser toutes les interactions du modèle. En l’occurrence, c’est l’utilisateur qui décide les interactions qu’il désire modéliser. Dans le cas présent, seule l’interaction entre le traitement et le sexe est modélisée.

• Le 1 de l’effet aléatoire indique que l’effet aléatoire est associé à l’intercept.

• La barre verticale sera la symbolique utilisée pour préciser que la variable est un facteur aléatoire.

4.3.1. La fonction ‘lmer’

La fonction ‘lmer’ du package ‘lme4’ a 4 arguments utiles qui sont décrits dans le Tableau 3.

Tableau 3. Description des arguments de la fonction ‘lmer’.

Argument	Description	Illustration
formula	il s’agit du modèle (effets fixes et aléatoires)	lmer(weight ~ Treatment + sex + Lsize +Treatment:sex + (1 | Litter),
data	Nom du jeu de données	data = rat,
REML A	lgorithme (argument logique TRUE/FALSE)	REML=T,
na.action	Manière de gérer les valeurs manquantes	na.action=na.omit,

De manière concrète, on utilise la fonction comme suit :

Modele.lmer<-lmer(  weight ~ Treatment + sex + Lsize +Treatment:sex+(1 |Litter), data = rat, 
             REML=T) 

Remarquez que l’effet aléatoire est entre parenthèse et directement intégré à l’équation du modèle.

4.4. Hypothèses

Dans notre exemple, on peut imaginer que les rats issus d’une même portée vont avoir un poids plus similaires que les rats issus de portées différentes. Dès lors nous sommes dans une situation de modèle multiniveau où un rat appartient à une portée. Donc, il y aurait deux sources de variance : La variance entre les portées qui sont au niveau 2, et la variance entre les rats issus d’une même portée qui sont au niveau 1.

4.4.1. Niveau 1

\[ Poids_ij = b_{0j} + b_{1j} × SEXE1_{ij} + \epsilon_{ij}\]

Le niveau 1 considère que le poids pour un rat i appartenant à une portée j suit un modèle typique d’une d’anova à une facteur avec un intercept \(b_{0j}\) qui est spécifiquement défini par la portée. De même, l’effet du sexe est spécifique à la portée j du rat. C’et ce que nous avons montré en comparant le modèle linéaire classique avec le modèle linéaire mixte.

Donc, les deux coefficients de régressions sont non observables (raison pour laquelle les coefficients sont en lettre romaines ici), et sont définis comme dépendant (i.e., étant fonction de) du niveau 2 du modèle.

4.4.2. Niveau 2

\[ b_{0j}= \beta_0+\beta_1 × TRAIT_{1j} + β2 × TRAIT_{2j} + β4 × Taille.Portée + u_j \] Et

\[ b_{1j} = \beta_3 + \beta5 × TRAIT1_{j} + \beta_6 × TRAIT2 \]

Ainsi, l’intercept pour l’effet du sexe (le niveau 1) sur le poids à la naissance dépend du traitement auquel a été soumis la mère ainsi que de la taille de la portée. De même, on peut faire l’hypothèse que l’effet spécifique du sexe dépend non seulement de l’effet spécifique du sexe (\(\beta_3\)) mais aussi du type de traitement auquel la mère a été soumise (on peut imaginer que le traitement n’a pas le même impact sur les males et sur les femelles).

5. Mise en application

5.1. Décrire les données

Avant de commencer à analyser les données, il est toujours utile de regarder les statistiques descriptives pour : - Identifier si a priori, il y a des différences entre les conditions ; - S’assurer que les données sont compatibles avec ce qu’elles sont censées être (en vérifier le minimum et le maximum)

Dans le package “psych”, le fonction ‘describeBy’ rempli cette fonction.

describeBy(x=rat$weight, group=list(rat$Treatment, rat$sex), mat=T )

##     item  group1 group2 vars  n     mean        sd median  trimmed
## X11    1 Control Female    1 54 6.116111 0.6851179  6.175 6.137045
## X12    2    High Female    1 32 5.851562 0.6001887  5.760 5.833846
## X13    3     Low Female    1 65 5.837538 0.4504964  5.840 5.823962
## X14    4 Control   Male    1 77 6.471039 0.7537880  6.410 6.450317
## X15    5    High   Male    1 33 5.918485 0.6909058  5.690 5.865926
## X16    6     Low   Male    1 61 6.025082 0.3803403  6.000 6.008980
##          mad  min  max range       skew    kurtosis         se
## X11 0.511497 3.68 7.57  3.89 -0.6327927  1.56777409 0.09323273
## X12 0.659757 4.48 7.68  3.20  0.4825407  1.05519051 0.10609938
## X13 0.385476 4.75 7.73  2.98  0.9281147  3.53662727 0.05587720
## X14 0.652344 4.57 8.33  3.76  0.1832662  0.01074036 0.08590211
## X15 0.637518 5.01 7.70  2.69  0.7237273 -0.42052370 0.12027125
## X16 0.370650 5.25 7.13  1.88  0.4353957  0.12413303 0.04869759

Il est par ailleurs toujours utile d’avoir une représentation visuelle d’un graphique.

data_summary <- function(x) {
  m <- mean(x)
  ymin <- m-sd(x)
  ymax <- m+sd(x)
  return(c(y=m,ymin=ymin,ymax=ymax))
}

# on crée le graphique en indiquant le jeu de données et les variables d'intérêt
p<- ggplot(data=rat, aes(x=Treatment, y=weight, fill=sex) )  
# on indique la forme du graphique, un graphe violon en l'occurrence 
p<-p+ geom_violin() 
 # ici on choisit le pallette de couleur (automatique pour ne pas devoir spécifier manuellement)
p<-p+scale_fill_brewer(palette="PRGn")
# on indique où la légende doit apparaître
p<-p + theme(legend.position="right") 
# on affiche chaque point individuel pour identifier la présence éventuelle de valeurs influentes 
p<- p+ geom_dotplot(binaxis = "y", stackdir = "center", position = "dodge",dotsize=1/4)
# on ajoute la moyenne par un point rouge et l'écart-type par un trait rouge. 
p<-p + stat_summary(fun.data=data_summary,geom="pointrange", color="red", size=0.50,position=position_dodge(0.9))
p

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Figure 3. Graphiques représentants les statistiques descriptives de la base de données ‘rat’

On remarque sur le graphique que certaines femelles ont un poids qui s’éloigne assez fortement de la moyenne de leur groupe, ce qui amène à suspecter la préence de valeurs influentes.

5.2. Valeurs influentes et conditions d’application

Pour identifier la présence de valeurs influentes et le respect des conditions d’application, il faut comprendre la notion de résidu. Le résidu est la différence entre la valeur observée et la valeur prédite par le modèle. Le résidu peut être conditionnel ou inconditionnel.

Il existe des résidus conditionnelles et inconditionnels.

Le résidu conditionnel est la différence entre la valeur observée est la valeur prédite conditionnelle, c’est-à-dire la valeur prédite en prenant en compte l’effet aléatoire. Le résidu inconditionnel est le résidu quand on ne prend pas en compte l’effet aléatoire

Prenons comme illustration les 5 premiers rats du jeu de données. Voici les poids qu’ils avaient à leur naissance :

rat$weight[1:5] # 5 premières observations du jeu de données

## [1] 6.60 7.40 7.15 7.24 7.10

Selon notre modèle, les valeurs prédites pour ces rats sont :

Modele.lme$fitted[1:5] # 5 premières valeurs prédites 

## [1] 6.785207 6.785207 6.785207 6.785207 6.785207

Le résidu est la différence entre ces valeurs observées et ces valeurs prédites.

rat$weight[1:5]-Modele.lme$fitted[1:5]

## [1] -0.1852067  0.6147933  0.3647933  0.4547933  0.3147933

On peut les obtenir de la manière suivante :

Modele.lme$residuals[1:5]

## [1] -0.1852067  0.6147933  0.3647933  0.4547933  0.3147933

Pour ces résidus, on parle de résidus bruts. Ceux-ci ne sont pas parfaitement adaptés pour tester les conditions d’application car les variances entre les groupes pourraient ne pas être parfaitement égales et les résidus conditionnels pourraient être corrélés.

Il est donc préférable de “studentiser” les résidus. Il s’agit de diviser le résidu par l’écart-type de chacun des groupes. Une alternative à cette méthode est de diviser les résidus par l’écart-type de la variable dépendante. On parle dans ce cas de résidus de Pearson.

Les résidus studentisés sont le plus souvent utilisés.

On obtient ces résidus normalisés avec la fonction “residuals”. L’argument ‘normalized’ permet d’avoir les résidus studentisés. Pour avoir les résidus de Pearson, il suffit de remplacer l’argument par “pearson”.

residus<-residuals(Modele.lme, type="normalized")
residus[1:10]

##           1           1           1           1           1           1 
## -0.67038237  0.88363473  0.39800439  0.57283131  0.30087832 -1.75819434 
##           1           1           1           1 
##  0.06777575  0.20375225  0.80230408 -0.47976003

5.2.1. Valeurs influentes

5.2.1.1. Fonction ‘lme’

La manière de penser les valeurs influentes peut être différente entre la fonction ‘nlme’ et la fonction ‘lmer’

Pour la fonction ‘nlme’, il n’existe pas de fonction qui a été développée pour tester spécifiquement l’effet de valeurs influentes contrairement à la fonction ‘lmer’. La logique peut donc être quelque peu différente entre les deux fonctions.

Dans le premier cas, on commence par récupérer les résidus du modèle (#1), on ajoute les résidus au jeu de données dans la variable ‘residus’ (#2) et on va tester si certains résidus sont particulièrement éloignés à l’aide du test de Grubbs qui teste la probabilité que le résidu le plus éloigné soit aussi éloigné au regarde de la taille de l’échantillon (#3).

residus<-residuals(Modele.lme, type="normalized") #1
rat$residus<-residus #2 
grubbs.test(rat$residus, type = 10, opposite = FALSE, two.sided = FALSE)

## 
##  Grubbs test for one outlier
## 
## data:  rat$residus
## G = 6.14055, U = 0.88217, p-value = 3.971e-08
## alternative hypothesis: lowest value -5.88670114350002 is an outlier

En l’occurrence, au moins un des rats considéré comme une valeur influente.

Il est possible de nettoyer le jeu de données des valeurs influentes. Pour éviter de commettre une erreur, il est toujours préférable de ne pas écraser le jeu de données mais de travailler sur un nouveau jeu de données. On va donc commencer par stocker le jeu de données ‘rat’ dans un jeu de données appelé ‘rat.clean’

rat.clean<-rat

Dans un second temps, on va créer une boucle qui permet de tester si la valeur de la probabilité du test de Grubbs est inférieure à 0.05 (i.e., il y a au moins une valeur influente) et supprimer l’observation qui a le résidu maximum en valeur absolue.

# On réalise un boucle du type: tant que la probabilité est inferieure a 0.05, on continue. 
data.frame()->valeur.influentes
while(grubbs.test(rat.clean$residus, type = 10, opposite = FALSE, two.sided = FALSE)$p.value <0.05)  { 
      max<-which.max(abs(rat.clean$residus)) #cherche la valeur maximale qu'on stocke dans l'objet max                            # récupère les observations considérées comme influentes et les stocke                                                        
      valeur.influentes<-rbind(valeur.influentes,rat.clean[max, ])
      rat.clean<-rat.clean[ -max, ] # supprime la valeur maximale de rat.clean
}  

Ceci nous permet de savoir qu’il y a deux observations considérées comme influentes,

nrow(rat)-nrow(rat.clean)

## [1] 2

Ce qui représente 0.62% de l’échantillon total.

100*(nrow(rat)-nrow(rat.clean)) /nrow(rat)

## [1] 0.621118

Les deux observations influentes sont :

valeur.influentes

## # A tibble: 2 x 7
##   id    weight sex    Litter Lsize Treatment residus
##   <fct>  <dbl> <fct>  <fct>  <dbl> <fct>       <dbl>
## 1 66      3.68 Female 6          9 Control     -5.89
## 2 227     4.75 Female 18        15 Low         -4.48

5.2.1.2. Fonction ‘lmer’

La même logique peut être appliquée avec la fonction ‘lmer’. Cependant, les résultats ne seront pas exactement les mêmes car les variances hétérogènes ne peuvent pas être modélisées avec la fonction ‘lmer’ alors qu’on l’a modélisée avec la fonction ‘lme’ (voir les points 4.3.1. et 4.3.2.). On voit en effet, qu’il y a 4 valeurs influentes en l’occurrence, mais les deux qui ont été identifiées précédemment le sont encore ici.

residus2<-residuals(Modele.lmer, type="pearson",scaled=TRUE) 
rat$residus2<-residus2 
grubbs.test(rat$residus2, type = 10, opposite = FALSE, two.sided = FALSE)

## 
##  Grubbs test for one outlier
## 
## data:  rat$residus2
## G = 7.79038, U = 0.81035, p-value = 3.575e-14
## alternative hypothesis: lowest value -7.47250745924872 is an outlier

rat.clean2<-rat
data.frame()->valeur.influentes
while(grubbs.test(rat.clean2$residus2, type = 10, opposite = FALSE, two.sided = FALSE)$p.value <0.05)  { 
      max<-which.max(abs(rat.clean2$residus2)) #cherche la valeur maximale qu'on stocke dans l'objet max                            # récupère les observations considérées comme influentes et les stocke                                                        
      valeur.influentes<-rbind(valeur.influentes,rat.clean2[max, ])
      rat.clean2<-rat.clean2[ -max, ] # supprime la valeur maximale de rat.clean
}
valeur.influentes

## # A tibble: 4 x 8
##   id    weight sex    Litter Lsize Treatment residus residus2
##   <fct>  <dbl> <fct>  <fct>  <dbl> <fct>       <dbl>    <dbl>
## 1 66      3.68 Female 6          9 Control     -5.89    -7.47
## 2 56      5.02 Female 5         13 Control     -3.01    -3.88
## 3 227     4.75 Female 18        15 Low         -4.48    -3.17
## 4 60      8.33 Male   6          9 Control      2.35     3.01

Une autre manière de tester la présence de valeurs influentes est de calculer les mesures d’influences, comme la distance de Cook. Cette manière de faire est possible uniquement avec la fonction ‘lmer’, bien que cette approche soit bien plus longue.

Ainsi, on peut tester l’effet de la suppression d’une portée sur l’estimation des paramètres. Il est important de constater que une des rattes qui était considérée comme influente provenait de la portée numéro 6.

cd<-case_delete(Modele.lmer, group="Litter", 
            type="both")
cd1<-diagnostics(cd) # on stocke dans l'objet cd1 les résultats des indicateurs diagnostics. 

On peut par exemple prendre la distance de Cook. Celle-ci permet d’identifier que les portées numéro 6 et 9 semblent avoir un gros impact sur l’estimation des paramètres.

dotplot_diag(x = COOKSD, data=cd1[["fixef_diag"]],cutoff="internal",name="cooks.distance",
             modify="dotplot",ylab = "Cooks Distance") # le cutoff interne est spécifié en fonction de name

Il est possible d’enlever une portée du jeu de données pour ne pas avoir d’observations qui influencent exagérément les résultats de la manière suivante :

rat2<-rat[which(rat$Litter!=6),] # supprime les rats issus de la portée numéro 6

Exercice 4 : faites en sorte de supprimer la portée numéro 9 en plus de la portée numéro 6

Il est possible d’obtenir d’autres indicateurs d’influence comme la différence d’ajustement (DDFITS) ou encore le ratio de covariance. De manière générale, on considère que la distance de Cook et ou la différence d’ajustement sont les meilleurs indicateurs. A titre informatif, voici les indicateurs pour la différence d’ajustement, et le ratio de covariance.

# DDFITS
dotplot_diag(x = MDFFITS, data=cd1[["fixef_diag"]],cutoff="internal",name="mdffits", 
             modify="dotplot",ylab = "MDFFITS")

# COVRATIO
dotplot_diag(x = COVRATIO, data=cd1[["fixef_diag"]],cutoff="internal",name="covratio", 
             modify="dotplot",ylab = "COVRATIO")

Dans les modèles linéaires mixtes, on peut chercher à identifier des valeurs influentes à la fois sur l’estimation des paramètres dans l’effet fixe, mais également dans l’estimation des variances des effets aléatoires. On obtient ces informations de la manière suivante :

dotplot_diag(x = sigma2, data=cd1[["varcomp_diag"]],cutoff="internal",name="mdffits", 
             modify="dotplot",ylab = "sigma2")

Une observation influente diminue considérablement l’estimation de la variance (elle la diminue de plus de “-0.20”). Il est possible d’identifier la portée à l’origine de cette diminution en demandant à R de préciser quelles sont les portées qui diminue la variance de plus de -0.20.

cd1[["varcomp_diag"]][which(cd1[["varcomp_diag"]]$sigma2<(-0.2)),] 

##   IDS     sigma2        D11
## 6   6 -0.2769763 0.07493979

Note importante : identifier des valeurs influentes n’implique pas nécessairement leur suppression. Cela dépend du contexte, de la légitimité des valeurs influentes, et des raisons à l’origine de cette influence. Néanmoins, de manière générale, laissez des valeurs influentes dans le jeu de données entraîne une estimation biaisée des paramètres. En revanche, il faut décider a priori des critères de suppression des valeurs influentes et s’y tenir.

Pour les besoins didactiques, nous comparerons le jeu de données complet avec le jeu de données nettoyées. Cela permettra d’illustrer l’impact de ces valeurs influentes sur plusieurs paramètres des modèles linéaires mixtes.

5.2.2. Vérifier les conditions d’application

Les MLM sont dépendants de deux conditions d’application : les résidus doivent être distribués en suivant une distribution normale et les variances aléatoires doivent également suivre une distribution normale.

Par ailleurs, le facteur aléatoire doit être indépendant du résidu.

Nous commencerons par tester la normalité, et ensuite l’indépendance entre le facteur aléatoire et le résidu.

5.2.2.1 La normalité

Tester la normalité du résidu est assez simple grâce au test de Shapiro-Wilk.

n1<-shapiro.test(rat$residus)
n2<-shapiro.test(rat.clean$residus)
n3<-shapiro.test(rat.clean2$residus)
r<-data.frame(W=c(n1$statistic, n2$statistic, n3$statistic),
              p=c(n1$p.value, n2$p.value, n3$p.value))
dimnames(r)[[1]]<-c("jeu de données complet", "jeu de données nettoyées lme", "jeu de données nettoyées lmer")
kable(pandoc.table(r, style='simple',split.tables=150))

	W	p
jeu de données complet	0.9465	2.052e-09
jeu de données nettoyées lme	0.9966	0.7417
jeu de données nettoyées lmer	0.9981	0.9763

## Warning in kable_markdown(x, padding = padding, ...): The table should have
## a header (column names)

|| || || ||

Il est également toujours utile de regarder la distribution des résidus ainsi que le QQplot. On obtient la distribution des résidus de la manière suivante :

p1<-ggplot(rat, aes(x=residus))+geom_histogram(aes(y=..density..))
p1<-p1+ stat_function(fun = dnorm, colour = "red",
                      args = list(mean = mean(rat$residus, na.rm = TRUE),
                                  sd = sd(rat$residus, na.rm = TRUE)))
p1<-p1+theme(plot.title = element_text(size = 12))+labs(x = "Distribution du résidu")
print(p1)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ce graphique s’interprète assez aisément : il s’agit de déterminer si les batonnets noirs s’ajustent à la distribution rouge. Dans notre exemple, on constate qu’il y a une asymétrie. Qu’en est-il lorsqu’on travaille sur les données nettoyées ?

Exercice 5 : Faites la représentation des résidus pour les données nettoyées avec la fonction lme et avec la fonction lmer

Le QQplot s’obtient de la manière suivante :

p2<-ggplot(rat, aes(sample=residus))+stat_qq() 
p2<-p2+theme(plot.title = element_text(size = 12))+ggtitle("QQplot")
print(p2)

Si les résidus de l’échantillon suivent une distribution normale théorique, les points doivent formés une ligne droite. On se rend compte que ce n’est pas le cas sur les extrémités.

Exercice 6 : Faites le QQplot pour les données nettoyées avec la fonction lme et avec la fonction lmer

En revanche, cette approche ne teste pas si la distribution de la variable aléatoire suit une distribution normale. Il est possible de tester la normalité de l’effet aléatoire par d’autres moyens mais uniquement lorsqu’on utilise le package lmer. La première étape consiste à modéliser le facteur aléatoire avec la fonction ‘lmer’. Pour indiquer à R que seule la constante doit être modélisée, on doit écrire le modèle de la manière suivante : \[ VD \sim 1 \]

Dans notre cas, on veut la constante et l’effet aléatoire. On les obtient donc de la manière suivante

aleatoires<-lmer(weight~1+(1|Litter), data=rat)

Grâce à la fonction ‘profile’, nous allons pouvoir dans un second temps faire un graphique de la vraisemblance de la variance.

pr01 <- profile(aleatoires)
xyplot(pr01, aspect = 1.3, layout=c(3,1))

Si la ligne est droite, cela signifie que la normalité est respectée, si elle est courbe, cela signifie que la normalité n’est pas respectée.

Dans notre cas, la variance du résidu est plutôt distribuée en suivant une distribution normale, mais la variance de la variable aléatoire ne suit pas totalement une distribution normale.

Ce phénomène sera plus manifeste lorsqu’on représente l’évolution de la vraisemblane en valeur absolue.

xyplot(pr01, aspect = 1.3, layout=c(3,1), absVal=T) 

## permet d'observer de manière plus efficace les intervalles de confiance

Sur le graphique du milieu et de droite, le V est symétrique. En revanche, sur celui de gauche, il ne l’est pas. Cela signifie que, si on estimait la variance un peu en-dessous de sa valeur, la vraisemblance diminuerait rapidement (sur le graphique, cela augmente, mais il faut garder à l’esprit que nous utilisons -2*log vraisemblance, donc les valeurs faibles signifie que la vraisemblance est élevée), alors que sur la partie droite, l’augmentation serait plus lente.

Exercice 7 : testez la normalité de l’effet aléatoire pour les données nettoyées par la fonction lmer.

5.2.2.2 Indépendance entre la variable aléatoire et le résidu

Le profil de la déviance fournit aussi des informations sur la sensibilité du modèle à des changements dans les paramètres. en effet, on peut représenter un graphique de coordonnées parallèles grâce à la fonction ‘splom’

splom(pr01)

Dans ce graphique, la diagonale représente les paramètres estimés. Ainsi, en haut à gauche, le graphique représente la manière dont la variance de l’effet aléatoire (sig01) va influencer la constante ; au milieu à gauche, cela indique l’influence entre l’effet aléatoire et le résidus ; et au haut au milieu, l’influence du résidu sur l’estimation de l’intercept.

Le graphique est symétrique, avec comme nuance que sur la partie inférieure, les résidus sont standardisés. Ainsi, il ne sont pas influencés par des problèmes d’échelle. Ce sont donc les trois graphiques sous la diagonale que nous observons.

Le contour des cercles représente les intervalles de confiance à un seuil fixé.

Les deux droites qui sont plus ou moins verticale et plus ou moins horizontale, représentent la manière dont un paramètre est estimé pour une valeur conditionnelle de l’autre. Dans notre exemple, le graphique au milieu à droite présente deux axes qui sont perpendiculaires. Cela signifie que les paramètres estimés ne s’influencent pas. En l’occurrence, la résidu n’est pas corrélé à l’intercept et leur estimations sont indépendantes.

Pour le graphique en bas à droite, on constate que, lorsqu’on se déplace sur l’axe horizontale c’est-à-dire que la valeur de la constante change, la valeur de la variance de l’effet aléatoire augmente. Les deux paramètres ne sont pas indépendants. Ceci est tout à fait normal puisque la variable aléatoire dépend de la constante, comme on le modélise explicitement dans la fonction ‘lme’.

Le dernier graphique (en bas au milieu) être plus important : comme les axes sont plutôt orthogonaux entre l’effet aléatoire et le résidu, cela indique que ces paramètres sont relativement indépendants l’un de l’autre.

5.3. Tester l’intérêt de l’effet aléatoire

Lorsque nous avons comparé le modèle linéaire au MLM, nous avons montré que, dans la pire des situations, le MLM n’était rien d’autre qu’un modèle linéaire. Il est donc raisonnable de commencer par se questionner sur l’intérêt de réaliser un MLM plutôt qu’un modèle linéaire. Pour cela, il faut tester l’intérêt de l’effet aléatoire.

Il existe plusieurs manières de s’assurer l’intérêt d’un effet aléatoire. Le premier est vérifier si la vraisemblance du modèle est améliorée lorsque l’effet aléatoire est présent en comparaison à un modèle où il n’est pas présent. Comme l’hypothèse porte sur effet aléatoire, l’algorithme qui sera utilisé sera le maximum de vraisemblance restreint.

Une seconde possibilité est de faire une analyse graphique de l’effet aléatoire.

Une troisième possibilité est de calculer le coefficient de corrélation intraclasse (ICC).

5.3.1. L’intérêt de l’effet aléatoire avec la fonction lme

On ne peut pas modéliser un MLM sans effet aléatoire. Il n’est donc pas possible de comparer directement avec la fonction lme un modèle avec le facteur aléatoire et un modèle où il ne serait pas présent.

Il est donc nécessaire de contourner cette difficulté. En réalité, il suffit de comparer les résltats du MLM avec les résultats d’un modèle des moindres carrés généralisés. la fonction gls dans R permet de réaliser un modèle des moindres carrés généralisés.

Cette fonction s’utilise exactement de la même manière que la fonction lme, excepté le fait qu’il n’y a pas de facteur aléatoire à spécifier.

Modele.GLS <- gls(weight ~ Treatment + sex + Lsize +Treatment:sex, data = rat, method = "REML",
                  weights = varIdent(form = ~1 | Treatment))

Une fois ce modèle réalisé, on peut le comparer avec le modèle lme que nous avons modélisé plus haut. Pour cela, on utilise la fonction anova en entrant les noms des modèles qu’on désire comparer les uns à côté des autres en les séparant d’une virgule.

anova(Modele.GLS, Modele.lme)

##            Model df      AIC      BIC    logLik   Test L.Ratio p-value
## Modele.GLS     1 10 489.9346 527.4604 -234.9673                       
## Modele.lme     2 11 381.8847 423.1630 -179.9423 1 vs 2  110.05  <.0001

Comme l’effet est significatif, cela suggère que l’effet de la portée mérite d’être gardée dans le modèle. Il faut néanmoins noter que la valeur de la probabilité qui est fournie ici n’est pas correcte parce qu’on teste l’hypothèse selon laquelle la variance de l’effet aléatoire de la portée est égale à 0. Cependant, une variance égale à 0 est une limite de l’espace de paramètre du calcul de la variance (i.e., une variance ne peut pas être égale à 0). Pour cette raison, le \(\chi^2\) est un mélange entre les 2 distributions avec un poids égal à 0.5. Il faudrait donc la diviser par deux. Néanmoins, divisée ou non par deux, cela ne change rien à l’interprétation.

Remarquez que deux modèles sont comparés en l’occurrence. En réalité, il est possible de comparer plusieurs modèles en même temps.

Note, on aurait obtenu les mêmes résultats si on avait précisé que la variable aléatoire était chacun des bébés rats qui sont nés (la variable ‘id’ dans notre modèle).

Modele.lme.id <- lme(weight ~ Treatment + sex + Lsize +Treatment:sex, random = ~ 1 | id,data = rat, 
                     method = "REML",weights = varIdent(form = ~1 | Treatment))
anova.out<-anova(Modele.lme, Modele.lme.id)
2*diff(anova.out$logLik)

## [1] -110.05

On obtient ce résultats en multipliant par deux la différence entre les deux logarithme de vraisemblance.

2*(-234.9673-(-179.9423))

## [1] -110.05

La raison pour laquelle le ratio de vraisemblance n’est pas calculé avec cette seconde procédure est que les deux modèles ne sont pas emboîtés.

5.3.2. L’intérêt de l’effet aléatoire avec la fonction lmer

Le package ‘lme4’ permet de directement tester l’intérêt de l’effet aléatoire grâce à la fonction ranova avec la probabilité calculée correctement.

ranova(Modele.lmer)

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## weight ~ Treatment + sex + Lsize + (1 | Litter) + Treatment:sex
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>          9 -200.55 419.10                         
## (1 | Litter)    8 -245.25 506.51 89.406  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exercice 8 : sans écraser les modèles stockés dans la mémoire de R, montrez que le ratio de vraisemblance pour l’effet aléatoire le même avec la fonction lmer et lme si on ne modélise pas les variances hétérogènes.

5.3.3. L’intérêt de l’effet aléatoire par représentation graphique

Ces résultats peuvent être aussi abordé sous un angle visuel en regardant si les effets aléatoires sont proches de 0 ou s’ils présentent une dispersion. Cette approche est cependant possible uniquement avec le package ‘lme4’.

Dans le graphique représenté à la Figure XX, on se rend compte que la variance à l’intérieur des portées (les barres horizontales) est plus faible que la variance entre les portées (manière dont les points bleus se dispersent)

dotplot(ranef(Modele.lmer, condVar=T))

## $Litter

Figure XX. Représentation de l’effet aléatoire : la dispersion entre les portées est plus grandes qu’à l’intérieur des portées.

5.3.4. L’intérêt de l’effet aléatoire calcul de l’ICC

Habituellement, l’ICC est une mesure décrivant l’homogénéité des réponses sur la variable dépendante à l’intérieur d’un cluster ou d’une unité d’analyse. Dans le contexte du modèle multiniveau, il s’agit d’une mesure permettant d’estimer l’homogénéité des réponses à l’intérieur d’un cluster (la portée dans notre cas). Ce coefficient de corrélation intraclasse sera élevé si la variance totale de l’effet aléatoire est nettement plus grande que la variance résiduelle.

Une autre manière de comprendre cet ICC est de le considérer comme la corrélation entre les paires d’individus à l’intérieur d’un cluster. La corrélation à l’intérieur de la portée doit être plus grande qu’entre différentes portées. Dans ce cas-ci, l’ICC peut être négatif. En général, on préfère la première définition.

Plus cet ICC sera élevé, et plus il sera pertinent de favoriser les modèles linéaires mixtes par rapport à l’ANOVA.

Pour calculer l’ICC, il faut commencer par modéliser un modèle où il n’y a que l’intercept et le facteur aléatoire. Nous allons le faire pour la fonction lme :

Modele.lme0<-lme(weight~1, random=~1|Litter, data=rat)

Comme le coefficient de corrélation intraclasse est défini par la formule :

\[ ICC_{portée}= \frac {\sigma^2_{portée}}{\sigma^2_{portée}+\sigma^2} \]

Cela signifie qu’il faut identifier ces composantes de variances dans le modèle que nous venons d’ajuster. Dans le package ‘nlme’, on obtient aisément ces informations grâce à la fonction VarCorr

VarCorr(Modele.lme0)

## Litter = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept) 0.3003700 0.5480602
## Residual    0.1963076 0.4430661

On obtient donc le coefficient de corrélation intraclasse par le calcul suivant :

icc.calcul<-0.3003700 /(0.3003700 + 0.1963076)
icc.calcul

## [1] 0.6047585

Avec un ICC à \(0.6047585\), on peut considérer que le poids des rats est bien plus homogène à l’intérieur des portées qu’entre les portées.

Exercice 9 : réalisez le modèle avec la fonction lmer où seul la constante et le facteur aléatoire sont modélisés. Stocker ce modèle dans un objet appelé ‘Modele.lmer0’. Identifiez dans la sortie de résultats les éléments de variance (ou écart-type) qui vont servir au calcul de l’ICC

En l’occurrence, nous avons calculé l’ICC à la main, mais il existe évidemment une fonction R qui permet de calculer l’ICC automatiquement, que ce soit pour la fonction lme ou pour la fonction lmer.

Dans le premier cas, on va utiliser la fonction ICC1.lme. Cette fonction a trois arguments qui sont dans l’ordre : la variable dépendante, le facteur aléatoire, et le jeu de données.

ICC1.lme(weight, Litter, rat)

## The following objects are masked _by_ .GlobalEnv:
## 
##     residus, residus2

## [1] 0.6047585

En ce qui concerne la fonction icc du package ‘sjstats’, elle calcule l’ICC pour n’importe quel modèle ajusté avec une fonction du package ‘lme4’.

sjstats::icc(Modele.lmer0)

## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.605
##   Conditional ICC: 0.605

Remarquez que la fonction icc est précédée par “sjstats::”. Cela signifie que la fonction icc qui doit être utilisée doit provenir du package ‘sjstats’. On utilise cette procédure quand il y a une possibilité de conflit (i.e., avoir le même nom) entre deux fonctions issues de packages différents.

5.5. Tester les effets fixes

Jusqu’à présent, nous n’avons pas encore testé l’apport des effets fixes dans le modèle. Il existe deux grandes approches possibles pour tester les effets fixes dans les modèles linéaires mixtes : réaliser une table d’anova et tester l’évolution du maximum de vraisemblance.

Une alternative à ces deux approches est d’utiliser les critères informations (i.e., AIC et BIC). Plus d’informations sur ces critères d’informations sont fournies dans le supplément à la formation.

5.5.1. L’amélioration du maximum de vraisemblance

Comme l’objectif est de tester des effets fixes, l’alogithme qui doit être utilisé est le maximum de vraisemblance (et non le maximum de vraisemblance restreint). L’utilisation du ratio de vraisemblance n’est possible qu’avec des modèles qui sont emboîtés.

Considérons les modèles A et B, qui sont deux modèles en compétition. Le modèle A va être emboîté dans le modèle B si c’est un “cas particulier” du modèle B, c’est-à-dire si le modèle A est un sous-espace d’un modèle B plus général. En d’autres termes,

Pour le formuler autrement, les paramètres du modèle emboîté peuvent être obtenus en imposant des contraintes sur les paramètres du modèle plus général (Imposer une contrainte serait par exemple que les paramètres valent 0, ou que les moyennes sur un facteur soient égales entre elles).

Donc, le modèle de référence “contient” les paramètres qui doivent être testés, alors que le modèle emboîté non.

A doit être le même modèle que B avec des contraintes supplémentaires (i.e., de nouveaux facteurs à tester).

Il faut donc ajuster plusieurs modèles et tester si l’ajustement est amélioré entre les modèles. Il existe deux logiques de construction des modèles dans les MLM : commencer par le modèle le plus complet et supprimer des facteurs ayant un moindre intérêt ou commencer par le modèle le plus épuré et ajouter progressivement des facteurs. Pour plus d’informations sur ces deux stratégies, vous pouvez vous reporter à la section correspondante du supplément à la formation.

Peu importe la stratégie utilisée, il est important d’avoir à l’esprit que l’ordre d’entrée des variables est particulièrement important. En effet, la vraisemblance du modèle de référence est testée par rapport au modèle emboîté. Par exemple, si on a un modèle initial qui est \(VD \sim VI_1\) et qu’on veut tester l’impact d’une seconde variable indépendante dans un modèle initial (\(VD \sim VI_1 + VI_2\)), l’ajoute de la seconde variable va être testé par rapport au premier modèle, et non par rapport à la variable dépendante uniquement. Ainsi, l’ordre des variables est extrêmement important, en particulier si on veut contrôler l’impact de certaines variables.

Dans notre exemple sur les rats, la taille de la portée et le sexe sont des variables de contrôle, le traitement est la variable dont on veut tester l’impact et l’interaction est un effet dont on peut imaginer que l’effet du traitement n’est pas le même en fonction du sexe du bébé rat. Il s’ensuit qu’en utilisant une stratégie ascendante, on doit modéliser 4 modèles que nous allons illustrer à l’aide de la fonction lmer : 1) Le modèle avec la constante uniquement.

Modele.lmer1<- lmer(weight ~  1 + (1 | Litter),data = rat, REML = F)

2. Le modèle avec les covariables (i.e., taille de la portée et sexe)

Modele.lmer2<- lmer(weight ~  sex + Lsize + (1 | Litter),data = rat, REML = F)

3. Le modèle 2 auquel on ajoute l’effet du traitement

Modele.lmer3 <- lmer(weight ~  sex + Lsize +Treatment + (1 | Litter),data = rat, REML = F)

4. Le modèle 3 auquel on ajoute l’effet de l’interaction

Modele.lmer4 <- lmer(weight ~  sex + Lsize +Treatment +Treatment:sex+ (1 | Litter),data = rat, REML = F)

Il est possible à présent de comparer l’ensemble de ces modèles à l’aide de la fonction anova

anova(Modele.lmer1,Modele.lmer2,Modele.lmer3,Modele.lmer4)

## Data: rat
## Models:
## Modele.lmer1: weight ~ 1 + (1 | Litter)
## Modele.lmer2: weight ~ sex + Lsize + (1 | Litter)
## Modele.lmer3: weight ~ sex + Lsize + Treatment + (1 | Litter)
## Modele.lmer4: weight ~ sex + Lsize + Treatment + Treatment:sex + (1 | Litter)
##              Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)
## Modele.lmer1  3 472.12 483.45 -233.06   466.12                          
## Modele.lmer2  5 407.48 426.35 -198.74   397.48 68.6434      2  1.242e-15
## Modele.lmer3  7 392.79 419.21 -189.39   378.79 18.6955      2  8.716e-05
## Modele.lmer4  9 395.81 429.78 -188.91   377.81  0.9721      2      0.615
##                 
## Modele.lmer1    
## Modele.lmer2 ***
## Modele.lmer3 ***
## Modele.lmer4    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On constate que les covariables et l’effet de traitement améliorent la vraisemblance du modèle, ce qui n’est pas le cas de l’interaction. Ainsi, on pourrait garder un modèle plus parcimonieux en supprimant l’interaction du modèle.

Exercice 10 : Appliquez la même logique de construction du modèle avec la fonction lme

5.5.2. La table de l’ANOVA

La table de l’ANOVA peut être utilisée pour les tests d’hypothèses linéaire relatives à de multiples effets fixes.

On obtient la valeur du F en utilisant à nouveau la fonction anova, mais en n’e précisant n’indiquant qu’un seul modèle à l’intérieur de la fonction. Pour la fonction lme du package ‘nlme’, par rapport à une anova classique, le type de F correspond au calcul de la somme de carrés de type I, c’est-à-dire séquentielle. En d’autres termes, l’effet testé est conditionnel à l’ensemble des effets testés auparavant dans le modèle. On obtient la table de l’anova de la manière suivante :

Modele.lme4 <- lme(weight ~  sex + Lsize +Treatment +Treatment:sex, random=~1 | Litter,data = rat, method = "ML")
anova(Modele.lme4)

##               numDF denDF   F-value p-value
## (Intercept)       1   292 10362.347  <.0001
## sex               1   292    53.234  <.0001
## Lsize             1    23    35.681  <.0001
## Treatment         2    23    13.349  0.0001
## sex:Treatment     2   292     0.476  0.6214

En l’occurrence, le calcul du F est une estimation, qui ne suit pas parfaitement la distribution F,ce qui implique que le déterminer les degrés de liberté de manière similaire à ce qui est fait dans le modèle linéaire classique (ce qui est fait dans ‘nlme’) est critiquée. En réalité, il est préférable d’estimer les degrés de liberté en prenant en compte la présence d’effets aléatoires et le fait que les résidus puissent être corrélés, ce qui n’est pas le cas avec la fonction lme, mais c’est le cas pour la fonction lmer. Pour être précis, les probabilités ne sont pas directement fournies par le package ‘lme4’ en raison de la difficulté dans l’estimation des degrés de liberté que nous venons d’aborder. En revanche, il est possible de les obtenir grâce au package ‘lmerTest’.

Ce package propose plusieurs méthodes pour l’estimation de ces degrés de liberté, dont la méthode de Satterthwaite qui est celle que nous allons utiliser.

anova(Modele.lmer4, type=1) # fonctionne parce que les packages ont été chargés auparavant

## Type I Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## sex           8.1705  8.1705     1 321.94 50.5228 7.632e-12 ***
## Lsize         4.5672  4.5672     1  35.01 28.2416 6.189e-06 ***
## Treatment     4.4267  2.2134     2  28.46 13.6865 6.846e-05 ***
## sex:Treatment 0.1575  0.0788     2 306.83  0.4871    0.6149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Remarquez que les degrés de liberté ont une virgule, cela est dû à l’estimation de Satterthwaite. Il ne s’agit pas d’une erreur. Par ailleurs, un autre avantage de ce package est qu’on peut choisir le type de somme des carrés qui est utilisée dans le calcul en modificant l’argument “type”. En effet, dans les résultats qui précède, on peut voir que la somme des carrés qui a été choisie est une somme des carrés de type 1.

Si nous ajustons notre modèle en choisissant un autre ordre dans les variables (en commençant par traitement par exemple), la valeur des F va changer. Ainsi, le F vaut 13.68 pour la variable “Treatment” dans le modèle ci-dessus, mais vaut 5.23 dans lorsque la variable ‘Treatment’ est entrée en premier comme dans le modèle ci-dessous.

Modele.lmer4b <- lmer(weight ~  Treatment +sex + Lsize +Treatment:sex+ (1 | Litter),data = rat, REML = F)
anova(Modele.lmer4b, type=1)

## Type I Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Treatment     1.6929  0.8465     2  27.753  5.2341   0.01178 *  
## sex           6.8761  6.8761     1 310.568 42.5189 2.834e-10 ***
## Lsize         8.7218  8.7218     1  38.269 53.9315 8.120e-09 ***
## Treatment:sex 0.1575  0.0788     2 306.835  0.4871   0.61490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En revanche, la valeur du F pour les variables n’est pas modifée par l’ordre d’entrée des variables lorsqu’on choisit une somme des carrés de type 3 et vaut 13.41.

anova(Modele.lmer4, type=3)[,5:6]

##               F value    Pr(>F)    
## sex           47.0447 3.810e-11 ***
## Lsize         53.3621 9.049e-09 ***
## Treatment     13.4163 7.797e-05 ***
## sex:Treatment  0.4871    0.6149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(Modele.lmer4b, type=3)[,5:6]

##               F value    Pr(>F)    
## Treatment     13.4163 7.797e-05 ***
## sex           47.0447 3.810e-11 ***
## Lsize         53.3621 9.049e-09 ***
## Treatment:sex  0.4871    0.6149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tout comme lorsque les effets ont été testés par l’amélioration de la vraisemblance, il apparaît que l’interaction n’apporte pas d’explication supplémentaire au modèle. Pour les prochains modèles, elle ne fera donc plus partie du modèle à tester.

5.5.3. Les outils d’aide au diagnotic

Une autre manière de déterminer quels facteurs doivent être maintenus ou supprimer est d’utiliser des outils d’aide au diagnostic. La fonction step du package ‘lmerTest’ permet de rendre ce service. Il s’agit de lui présenter le modèle saturé (i.e., avec tous les effets) et il va déterminer quels sont les facteurs aléatoires et les effets fixes qu’il faudrait supprimer. En l’occurrence, il n’y a pas d’effet fixe qu’il faudrait supprimer. En revanche, le chiffre 1 dans la colonne “Eiminated” indique que l’interaction n’est pas indispensable au modèle et mérite qu’elle soit supprimée, contrairement aux autres effets.

st<-step(Modele.lmer4)
st

## Backward reduced random-effect table:
## 
##              Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                     9 -188.91 395.81                         
## (1 | Litter)          0    8 -231.14 478.27 84.461  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##               Eliminated Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)
## sex:Treatment          1 0.1575  0.0788     2 307.131  0.4871    0.6149
## sex                    0 9.2302  9.2302     1 304.361 56.9340 5.254e-13
## Lsize                  0 8.6775  8.6775     1  38.225 53.5249 8.909e-09
## Treatment              0 4.3844  2.1922     2  28.502 13.5221 7.419e-05
##                  
## sex:Treatment    
## sex           ***
## Lsize         ***
## Treatment     ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## weight ~ sex + Lsize + Treatment + (1 | Litter)

Cette fonction permet aussi d’obtenir une représentation graphique où ma couleur des batonnets indique le niveau de significativité.

plot(st) 

5.6. Les tailles d’effet

A présent que nous connaissans la significativité des différentes, il est souhaitable de connaître la taille de l’effet de notre modèle. Les normes APA encouragent à reporter de manière systématique la taille de l’effet dans les articles. Dans le cas du modèle linéaire mixte, deux stratégies peuvent être utilisées. Elles s’appuient toute deux sur le packages ‘MuMin’.

La première stratégie consiste à calculer les composantes de variances sur les effets fixes, aléatoires et le résidu afin de calculer le pourcentage de variance expliquée. Cette approche a du sens tant qu’on fait l’hypothèse d’une distribution normale (ce qui ne serait pas le cas par exemple sur une variable dépendante dichotomique où nous modéliserions la variable comme suivant une distribution binomiale).

Modele.final <- lme(weight ~  sex + Lsize +Treatment , random=~1 | Litter,data = rat, method = "ML")
r.squaredGLMM(Modele.final)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.

##            R2m       R2c
## [1,] 0.4346107 0.6238045

La première valeur est le \(R^2\) marginal et la seconde valeur est le \(R^2\) conditionnel. Le \(R^2\) marginal est le pourcentage de variance expliquée par les effets fixes et le \(R^2\) marginal est le pourcentage de variance expliquée par les effets fixes et aléatoires.

La seconde stratégie plus généraliste consiste à calculer un pseudo-R² en calculant l’évolution du maximum de vraisemblance. Les personnes familières avec la régression logistique connaisse le pseudo-R² de Nagelkerke. L’idée générale est de prendre la différence entre les deux maximum de vraisemblance et de déterminer le pourcentage d’évolution en divisant par le modèle nul.

 r.squaredLR(Modele.final)

## [1] 0.5456628
## attr(,"adj.r.squared")
## [1] 0.6345789

Cette fonctions renvoie deux valeurs. La première est le pourcentage d’amélioration du maximum de vraisemblance et la valeur du \(R^2_ajusté\) renvoie cette amélioration en appliquant la correction de Nagelkerke. Cette seconde valeur est censée être plus précise.

Exercice 11 : testez les R² pour les modèles ajustés par la fonction lmer.

5.7. Les contrastes

Quand on réalise une anova, on désire régulièrement pouvoir faire les contrastes de notre choix. En l’occurrence, les contrastes qui sont par défaut dans R, ce sont les contrastes qui permettent de faire une comparaison à une ligne de base.

En réalité, il est possible de modifier ces contrastes de différentes manières : en utilisant des matrices de contrastes par défaut, en créant la matrice de contastes, ou en testant des contrastes particulier.

5.7.1. Les matrices par défaut

Dans les matrices qu’on peut utiliser dans R par défaut, il y a 3 types de matrices :

1. Les comparaison à une ligne de base : il s’agit du contraste par défaut, ce qui est appelé “treatment” dans R.

Ainsi, si nous avons 3 modalités, on obtient cette matrice de contraste grâce à la fonction contr.treatment. Le premier argument est le nombre de modalité. L’argument “base” permet de spécifier la modalité qui représente la ligne de base. Dans cet exemple, il s’agit de la seconde modalité.

 contr.treatment(3,base =2, contrasts = TRUE, sparse = FALSE)

##   1 3
## 1 1 0
## 2 0 0
## 3 0 1

2. Les comparaison de type Helmert : il s’agit des contrastes orthogonaux.

Ainsi, si nous avons 3 modalités, on obtient cette matrice de contraste grâce à la fonction contr.helmert.

 contr.helmert(3)

##   [,1] [,2]
## 1   -1   -1
## 2    1   -1
## 3    0    2

Il est possible d’inverser cette matrice pour faire en sorte que la première modalité ait le coefficient 2 :

apply(contr.helmert(3), 2, rev)

##   [,1] [,2]
## 3    0    2
## 2    1   -1
## 1   -1   -1

3. Les comparaisons polynomiales. Lorsque notre hypothèse est que l’évolution entre les modalités n’est pas linéaire mais quadratique ou cubique, nous devons utiliser des contrastes polynomiaux qu’on obtient de la manière suivante :

 contr.poly(3)

##                 .L         .Q
## [1,] -7.071068e-01  0.4082483
## [2,] -7.850462e-17 -0.8164966
## [3,]  7.071068e-01  0.4082483

A présent, que les matrices disponibles par défaut dans R ont été définies, il faut attribuer une matrice de contrastes à un facteur. Imaginons qu’on veuille faire des contrastes du type Helmert sur la variable traitement, nous commençons par associer la variable ‘traitement’ à des contrastes de type Helmert :

 contrasts(rat$Treatment)<-contr.helmert(3) # le 3 représente le nombre de modalités du facteur

Nous réalisons notre modèle :

model.helmert <- lme(weight ~ Treatment + sex + Lsize , random = ~ 1 |Litter ,data = rat, method = "ML") 
# notez que l'interaction a disparu de la formule
summary(model.helmert)$tTable

##                     Value  Std.Error  DF     t-value       p-value
## (Intercept)  7.5201081812 0.22255344 294 33.79012353 2.966506e-103
## Treatment1  -0.4305150262 0.08512328  23 -5.05754753  4.044570e-05
## Treatment2   0.0004871619 0.04331222  23  0.01124768  9.911228e-01
## sexMale      0.3571708037 0.04770771 294  7.48664767  8.274578e-13
## Lsize       -0.1288262615 0.01774694  23 -7.25907099  2.179526e-07

Et nous pouvons obtenir la table des contrastes avec la fonction summary.

summary(model.helmert)$tTable

##                     Value  Std.Error  DF     t-value       p-value
## (Intercept)  7.5201081812 0.22255344 294 33.79012353 2.966506e-103
## Treatment1  -0.4305150262 0.08512328  23 -5.05754753  4.044570e-05
## Treatment2   0.0004871619 0.04331222  23  0.01124768  9.911228e-01
## sexMale      0.3571708037 0.04770771 294  7.48664767  8.274578e-13
## Lsize       -0.1288262615 0.01774694  23 -7.25907099  2.179526e-07

Dans cette table, on constate qu’il y a 2 contrastes liés au au traitement, un contraste lié au sexe et un contraste lié à la taille de la portée. Pour chaque contrastes, il y a les colonnes ‘value’, ‘Std.error’, ‘DF’, ‘t-value’ et ‘p-value’. Voici à quoi elles correspondent : - value : coefficient de la droite de régression. Quand il s’agit d’une covariable numérique, c’est la pente, quand il s’agit d’une variable catégorielle, c’est la différence de moyennes entre les modalités impliquées dans le contraste. Pour la constante, il s’agit de la moyenne de la variable dépendante. - Std.Error : il s’agit de l’erreur standard du coefficient de régression. - DF : les degrés de liberté pour l’effet testé. - t-value : le rapport entre value et Std.Error. Cette valeur permet de calculer la significativité du contraste au regard des ddl. - p-value : probabilité associée au contraste permettant de déterminer s’il existe ou non une différence significative.

En l’occurrence, tous les effets sont significatifs sauf le second contrastes du traitement.

Il est concevable que les numéros de contrastes ne soient pas bien clairs pour déterminer ce qui est comparé. Alors bien sûr pour le sexe où il n’y a que 2 modalités, les choses sont simples puisque le contraste compare nécessairement les males aux femelles. Cependant, pour le traitement, les choses sont moins claires. Un moyen de rendre les choses un peu claires est de reprendre la matrice de contrastes dans la variable ‘Treatment’ :

 contrasts(rat$Treatment)

##         [,1] [,2]
## Control   -1   -1
## High       1   -1
## Low        0    2

Ainsi, le premier contraste compare le groupe de contrôle au groupe ‘high’, tandis que le second contraste compare le groupe de ‘Low’ aux deux autres groupes. Ces contrastes ne sont pas les plus judicieux. Il aurait été préférable de comparer les deux groupes ayant eu un traitement avec le groupe de contrôle et ensuite, les deux groupes ayant eu un traitement entre eux.

Exercice 12 : en utilisant la fonction lmer, ajustez un modèle qui permet de comparer dans les contrastes le groupe de contrôle aux deux autres groupes et ensuite les groupes ayant eu un traitement entre eux. Maintenez dans le modèle l’interaction et identifiez la manière dont les contrastes d’interaction apparaîssent dans la table des t. Notez également le fait que les degrés de liberté ont été ajustés et ont donc une valeur avec virgule.

5.7.2. Les matrices définies par l’utilisateur

Avant de pouvoir spécifier soi-même les contrastes, il est nécessaire de bien comprendre la distinction entre une matrice et un vecteur, ainsi que la manière dont on crée des vecteurs et des matrices dans R.

5.7.2.1. Les matrices et les vecteurs

Un vecteur est un ensemble de données à une dimension (une ligne ou une colonne). Par exemple, si on peut stocker les chiffres de 1 à 10 dans l’objet ‘a’, Cette listes de chiffres est un vecteur. Pour créer un vecteur dans R, on utilise la fonction c.

a<-c(1:10)
a

##  [1]  1  2  3  4  5  6  7  8  9 10

Une matrice est un tableau de données dans lequel toutes les données sont de même nature. On dit que la matrice est carrée si le nombre de lignes correspond au nombre de colonnes, autrement, on dit qu’elle est rectangulaire. Pour créer une matrice dans R, on utilise la fonction matrix.

matrix(a,5) # une matrice rectangulaire

##      [,1] [,2]
## [1,]    1    6
## [2,]    2    7
## [3,]    3    8
## [4,]    4    9
## [5,]    5   10

a<-matrix(1:9,3) # une matrice carrée
a

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

On peut accéder à n’importe quel élément d’une matrice grâce à l’adressage matriciel. L’adressage matriciel consiste à donner les coordonnées entre crochets (‘[’ et ’]’) de l’information désirée. Dans l’adressage matriciel, il faut donner deux informations, séparées par une virgule. La première information est le numéro de la ligne et la seconde est le numéro de la colonne.

Par exemple, dans la table ‘a’ que nous venons de créer, on peut obtenir l’information dans la 2e ligne et dans la 2e colonne (i.e., 5) de la manière suivante :

a[2,2] 

## [1] 5

On peut obtenir aussi l’ensemble des informations d’une ligne ou d’une colonne en laissant un des arguments vides. Par exemple, on peut obtenir la 2e colonne d’une matrice :

a[,2] # fournit toute la seconde colonne

## [1] 4 5 6

Et on peut obtenir les informations dans la 3e lignes de la manière suivante :

a[3,] # fournit toute la 3e ligne

## [1] 3 6 9

5.7.2.2. Spécification de la matrice de contrastes par l’utilisateur

Comme les coefficients de contrastes sont des matrices, on peut faire notre matrice de manière manuelle. Chaque colonne doit représenter un vecteur de coefficient. Pour faciliter les choses, on va créer nos vecteurs séparément.

Imaginons que nous désirions appliquer des coefficients de type Helmert avec le 2 sur la second modalité. Nous pourrions faire notre matrice de la manière suivante :

cont1<-c(1,-2,1)
cont2<-c(1,0,-1)
cont1

## [1]  1 -2  1

cont2

## [1]  1  0 -1

On peut ensuite combiner ces vecteurs dans une matrice à l’aide de la fonction cbind qui va concatener une liste de vecteurs de sorte à ce que chacun représente une colonne.

cont<-cbind(cont1, cont2)
cont

##      cont1 cont2
## [1,]     1     1
## [2,]    -2     0
## [3,]     1    -1

La suite se réalise comme plus haut : on attribue ces contrastes à notre facteur. Ensuite, on ajuste le modèle et on obtient les coefficients par la fonction summary.

 contrasts(rat$Treatment)<-cont
Modele.manuel <- lme(weight ~ Treatment + sex + Lsize , random = ~ 1 |Litter ,data = rat, method = "ML") 
summary(Modele.manuel)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rat 
##        AIC      BIC    logLik
##   392.7857 419.2076 -189.3929
## 
## Random effects:
##  Formula: ~1 | Litter
##         (Intercept)  Residual
## StdDev:   0.2855393 0.4026421
## 
## Fixed effects: weight ~ Treatment + sex + Lsize 
##                    Value  Std.Error  DF  t-value p-value
## (Intercept)     7.520108 0.22255344 294 33.79012  0.0000
## Treatmentcont1  0.215501 0.05147587  23  4.18645  0.0004
## Treatmentcont2  0.214527 0.07017494  23  3.05703  0.0056
## sexMale         0.357171 0.04770771 294  7.48665  0.0000
## Lsize          -0.128826 0.01774694  23 -7.25907  0.0000
##  Correlation: 
##                (Intr) Trtmn1 Trtmn2 sexMal
## Treatmentcont1  0.352                     
## Treatmentcont2  0.035  0.010              
## sexMale        -0.067  0.009 -0.042       
## Lsize          -0.954 -0.418 -0.033 -0.045
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -7.5970401770 -0.4823977270  0.0009856047  0.5649441016  3.0646087033 
## 
## Number of Observations: 322
## Number of Groups: 27

IMPORTANT : les matrice de contrastes doivent impérativement respecter les règles de constructions des contrastes sinon, un message d’erreur sera renvoyé par R.

5.7.3. Les autres possibilités de contrastes.

Dans des modèles complexes, il arrive parfois qu’on désire obtenir les effets simples ou que tous les contrastes ne nous intéressent pas. Il arrive aussi qu’on veuille réaliser des contrastes d’interaction 2 à 2.

Dans tous ces cas, deux packages sont faits pour vous : le package ‘phia’ et le package ‘emmeans’. Le premier calcule les contrastes en utilisant le ratio de vraisemblance alors que le second calcule le t classique.

5.7.3.1. Le package ‘phia’

L’utilisation la plus simpliste de ce package pour tester les interactions est d’utiliser la fonction testInteractions en spécifiant uniquement le modèle. Pour avoir une interaction, nous utiliserons le modèle avec interaction, même si ce n’est pas le modèle final retenu. Dans ce cas, il va faire toutes les comparaisons 2 à 2.

testInteractions(Modele.lme4)

## Chisq Test: 
## P-value adjustment method: holm
##                                Value Df  Chisq Pr(>Chisq)
## Female-Male : Control-High -0.110009  1 0.7070          1
## Female-Male :  Control-Low -0.084144  1 0.6419          1
## Female-Male :     High-Low  0.025865  1 0.0381          1

A présent que nous avons les contrastes, nous pourrions vouloir les moyennes et erreurs-types ajustées, ainsi que leur représentation graphique afin de se faire une meilleure représentation de l’effet de chacune des variables

means<-interactionMeans(Modele.lme4)
means

##      sex Treatment adjusted mean std. error
## 1 Female   Control      6.201623  0.1063312
## 2   Male   Control      6.612176  0.1018926
## 3 Female      High      5.401904  0.1473636
## 4   Male      High      5.702448  0.1468528
## 5 Female       Low      5.818192  0.1044755
## 6   Male       Low      6.144601  0.1065743

plot(means, abbrev.levels=TRUE)

Pour éviter de multiplier l’erreur de première espèce (ou pour ne pas corriger la probabilité), on peut préciser la correction de la probabilité que nous désirons. Nous allons comparer les résultats sans correction de la probabilité et ceux avec une correction de holm :

testInteractions(Modele.lme4, adjustment="none") # sans correction

## Chisq Test: 
## P-value adjustment method: none
##                                Value Df  Chisq Pr(>Chisq)
## Female-Male : Control-High -0.110009  1 0.7070     0.4004
## Female-Male :  Control-Low -0.084144  1 0.6419     0.4230
## Female-Male :     High-Low  0.025865  1 0.0381     0.8452

testInteractions(Modele.lme4, adjustment="holm") # avec correction de Holm

## Chisq Test: 
## P-value adjustment method: holm
##                                Value Df  Chisq Pr(>Chisq)
## Female-Male : Control-High -0.110009  1 0.7070          1
## Female-Male :  Control-Low -0.084144  1 0.6419          1
## Female-Male :     High-Low  0.025865  1 0.0381          1

Face à des interactions, il arrive d’avoir besoin de tester les effets simples. Par exemple, nous pourrions vouloir savoir si l’effet du sexe se différencie sur chaque modalité de la variable ‘Treatment’. Pour cela, deux nouveaux paramètres doivent être ajoutés : la variable sur laquelle on veut les comparaisons (i.e., sexe) et celle sur laquelle nous voulons décomposer l’effet du sexe (i.e., traitement) :

testInteractions(Modele.lme4, adjustment="none", 
                 pairwise="sex", fixed="Treatment")

## Chisq Test: 
## P-value adjustment method: none
##                          Value Df   Chisq Pr(>Chisq)    
## Female-Male : Control -0.41055  1 31.8638  1.654e-08 ***
## Female-Male :    High -0.30054  1  7.6376   0.005716 ** 
## Female-Male :     Low -0.32641  1 18.5612  1.645e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ainsi,nous avons la différence entre les males et les femelles pour chacun des groupes pris séparément et tous ces contrastes sont significatifs.

La fonction testInteractions permet de tester d’autres contrastes que les contrastes 2 à 2. Dans ce cas, il faut préciser la matrice de contrastes désirée. Nous allons réutiliser notre matrice de contrastes utilisée plus haut sur traitement. Comme il n’est pas possible de tester en même temps des contrastes en comparaison 2 à 2 ainsi qu’une matrice déterminée, il faut supprimer l’argument ‘pairwise’ et ajouter l’argument ‘custom’. Pour que la fonction sache à quel facteur cette matrice de contraste doit être appliquée, il faut que la matrice ait le nom de la variable.

On doit donc commencer par créer une liste, il s’agit d’un format particulier de stockage des données dans R. Cette liste sera appelée “cont1”. Ensuite, on stocke la matrice de contrastes dans la liste avec comme nom la variable à laquelle s’applique les contrastes. En l’occurrence, la matrice ‘cont’ doit être appliquée à la variable ‘Treatment’.

cont1<-list()
cont1$Treatment<-cont
cont1

## $Treatment
##      cont1 cont2
## [1,]     1     1
## [2,]    -2     0
## [3,]     1    -1

testInteractions(Modele.lme4, adjustment="none", custom=cont1, pairwise="sex")

## Chisq Test: 
## P-value adjustment method: none
##                         Value Df  Chisq Pr(>Chisq)
## Female-Male : cont1 -0.135874  1 0.3167     0.5736
## Female-Male : cont2 -0.084144  1 0.6419     0.4230

Les résultats de la fonction ci-dessus indique que le \(\chi^2\) du contraste permettant de comparer la différence entre les males aux femelles sur l’effet du traitement entre ‘low’ et les deux autres groupes vaut 0.3167 et est non significatif. De même la différence entre les males et les femelles n’est pas plus marquée pour le groupe de contrôle par rapport au groupe ayant reçu un traitement élevé puisque le \(\chi^2\) vaut 0.6419 associé à une probabilité de 0.4230.

En réalité, il est possible de réaliser n’importe quel contraste avec cette fonction. Elle doit donc être utilisée avec prudence parce qu’on pourrait faire des contrastes tel que 5-2+3 et cela fonctionnerait bien que cela n’ait pas de sens.

5.7.3.1. Le package ‘emmeans’

Le second package avec lequel on peut réaliser des contrastes d’interaction, mesurer des effets simple et bien d’autres choses est le package ‘emmeans’.

Pour utiliser ce package, il faut travailler en deux étapes. La première consiste créer un objet ayant les caractéristiques demandées par ‘emmeans’. Cela se réalise en utilisant deux arguments de la fonction emmeans. Le premier argument est le modèle sur lequel on doit calculer les contrastes. Le package est compatible avec un grand nombre de packages et de modèle. Nous allons l’illustrer avec un modèle ajusté par le package ‘lmer’. Le second argument consiste à préciser les variables sur lesquelles on désire réaliser des contrastes. En l’occurrence, les deux variables d’intérêt sont ‘Treatment’ et ‘sex’. On va donc stocker ces informations dans un objet appelé “emmeans.out”.Cet objet fournit la liste des moyennes ajustées par combinaison de modalités, l’erreur-type ainsi que l’intervalle de confiance relatif à l’estimation de la moyenne.

Nous avons choisi plus haut d’appliquer une correction de Satterthwaite aux degrés de liberté pour le modèle lmer. Il faut donc commencer par préciser cette option dans le package “emmeans”. La fonction emm_options permet de préciser la manière dont les degrés de liberté seront calculés sur le modèle lmer. La fonction get_emm_option permet de s’assurer que les paramètres que nous avons choisi ont bien été pris en compte et on termine en utilisant la fonction emmeans

emm_options(lmer.df = "satterthwaite")
get_emm_option("lmer.df")

## [1] "satterthwaite"

emmeans.out<- emmeans(Modele.lmer4, ~ Treatment*sex, weights="show.levels")
emmeans.out

##  Treatment sex    emmean    SE   df lower.CL upper.CL
##  Control   Female   6.20 0.106 37.5     5.99     6.42
##  High      Female   5.40 0.147 39.3     5.10     5.70
##  Low       Female   5.82 0.104 33.9     5.61     6.03
##  Control   Male     6.61 0.102 31.7     6.40     6.82
##  High      Male     5.70 0.147 38.3     5.41     6.00
##  Low       Male     6.14 0.107 35.9     5.93     6.36
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95

Pour les personnes qui préfèrent avoir une représentation graphique pour se faire une idée de la manière dont les moyennes se répartissent, on l’obtenir aisément de la manière suivante :

plot(emmeans.out)

A partir de cet objet, on peut faire toutes les comparaisons 2 à 2 :

pairs(emmeans.out, adjust="holm") # toutes les comparaisons 2 à 2

##  contrast                      estimate     SE    df t.ratio p.value
##  Control,Female - High,Female     0.800 0.1821  40.3  4.392  0.0008 
##  Control,Female - Low,Female      0.383 0.1491  35.8  2.572  0.0721 
##  Control,Female - Control,Male   -0.411 0.0727 299.1 -5.645  <.0001 
##  Control,Female - High,Male       0.499 0.1816  39.5  2.748  0.0539 
##  Control,Female - Low,Male        0.057 0.1506  36.7  0.379  1.0000 
##  High,Female - Low,Female        -0.416 0.1785  39.1 -2.332  0.0759 
##  High,Female - Control,Male      -1.210 0.1795  38.1 -6.741  <.0001 
##  High,Female - High,Male         -0.301 0.1088 312.4 -2.764  0.0424 
##  High,Female - Low,Male          -0.743 0.1812  40.8 -4.098  0.0017 
##  Low,Female - Control,Male       -0.794 0.1460  32.9 -5.439  0.0001 
##  Low,Female - High,Male           0.116 0.1782  38.5  0.650  1.0000 
##  Low,Female - Low,Male           -0.326 0.0758 306.1 -4.308  0.0002 
##  Control,Male - High,Male         0.910 0.1791  37.4  5.079  0.0001 
##  Control,Male - Low,Male          0.468 0.1475  33.8  3.171  0.0258 
##  High,Male - Low,Male            -0.442 0.1808  40.0 -2.445  0.0759 
## 
## P value adjustment: holm method for 15 tests

On peut également faire les comparaisons 2 à 2 pour la variable ‘Treatment’ sur chacune des modalités de la variable ‘sex’, c’est-à-dire calculer les effets simples. Pour cela, on précise l’argument “by” qui sera la variable sur laquelle on va décomposer les effets simples.

pairs(emmeans.out, by = c("sex"), adjust="holm")

## sex = Female:
##  contrast       estimate    SE   df t.ratio p.value
##  Control - High    0.800 0.182 40.3  4.392  0.0002 
##  Control - Low     0.383 0.149 35.8  2.572  0.0288 
##  High - Low       -0.416 0.179 39.1 -2.332  0.0288 
## 
## sex = Male:
##  contrast       estimate    SE   df t.ratio p.value
##  Control - High    0.910 0.179 37.4  5.079  <.0001 
##  Control - Low     0.468 0.147 33.8  3.171  0.0064 
##  High - Low       -0.442 0.181 40.0 -2.445  0.0190 
## 
## P value adjustment: holm method for 3 tests

Pour terminer, on peut également réaliser d’autres matrices de contrastes grâce à la fonction contrast. Il est possible d’utiliser une grande variété de contrastes qui sont prédéfinis (et documentés dans le pdf du package ‘emmeans’ à la fonction contrast-methods) ou créer ses propres contrastes.

Nous allons commencer par illustrer un contraste prédéfini. Il semble y avoir du sens à comparer les groupes traités au groupe de contrôle. Pour cela, on peut utiliser une matrice prédéfinie dans ‘emmeans’ qui correspond à “trt.vs.ctrl1”. En réalité, les coefficients de ces contrastes sont obtenus par une fonction qui s’appelle trt.vs.ctrl1.emmc et qui a comme argument un argument qui s’appelle “ref”. Cet argument consiste à spécifier le numéro de la modalité qui est la modalité de référence. Dans le cas de l’exemple sur les rats, les modalités sont dans l’ordre “Control”,“High” et “Low”. La modalité de référence est donc la modalité 1.

Pour la variable ‘sex’, il n’y a qu’une comparaison possible. Le plus simple est donc d’utiliser une comparaison 2 à 2.

cont.int<-contrast(emmeans.out, interaction = c("trt.vs.ctrl1", "pairwise"), ref=1 )
cont.int

##  Treatment_trt.vs.ctrl1 sex_pairwise  estimate    SE  df t.ratio p.value
##  High - Control         Female - Male   0.1100 0.131 309 0.841   0.4011 
##  Low - Control          Female - Male   0.0841 0.105 303 0.801   0.4236

Dans notre cas, aucun des deux contrastes d’interaction n’est significatif (ce qui était attendu vu que l’interaction n’est pas significative).

Le package ‘emmeans’ permet aussi à l’utilisateur de spécifier manuellement les contrastes. Pour cela, il est nécessaire de comprendre la manière de rentrer les contrastes. Reprenons les contrastes d’interaction que nous venons de réaliser et intéressons-nous aux coefficients des contrastes à l’aide de la fonction coef. On constate que les deux premières colonnes correspondent à la combinaison des modalités des deux variables et les 2 dernières colonnes aux coefficients associés à chacune des combinaisons de modalités.

mes.coef<-coef(cont.int)
mes.coef

##   Treatment    sex c.1 c.2
## 1   Control Female  -1  -1
## 2      High Female   1   0
## 3       Low Female   0   1
## 4   Control   Male   1   1
## 5      High   Male  -1   0
## 6       Low   Male   0  -1

La fonction contrast peut directement utiliser des coefficients de contrastes. Cependant, elle requière uniquement des valeurs numériques. Pour le moment, ces coefficients sont dans les colonnes 3 et 4. Il faut donc préciser les colonnes dans lesquelles contrast va trouver ces coefficients. On peut le faire aisément dans R grâce à l’adressage matriciel.

contrast(emmeans.out,mes.coef[,3:4] )

##  contrast estimate    SE  df t.ratio p.value
##  c.1        0.1100 0.131 309 0.841   0.4011 
##  c.2        0.0841 0.105 303 0.801   0.4236

A présent que la manière dont contrast demande qu’on présente les coefficients, nous pouvons créer des vecteurs de contrastes. Quand il s’agit de tester les effets principaux, les choses sont encore relativement simples, mais elles peuvent se compliquer pour des contrastes d’interactions, a fortiori lorsqu’il y a plus que 2 facteurs. Une manière élégante et pas trop compliquée (tout en évitant de se tromper) de créer la bonne matrice de contraste consiste à commencer par créer les vecteurs pour chaque variable séparément. Pour l’illustrer, on va reprendre les contrastes que nous avions définis plus haut où on compare d’abord la modalité “High” aux deux autres et ensuite la modalité “Control” à la modalité “Low”. Pour améliorer la lisibilité, les noms qui sont données à ces deux vecteurs seront ‘Tr1’ et ‘Tr2’ pour ‘Traitement 1’ et ‘Traitement 2’. Pour la variable sexe, les coefficients valent nécessairement 1 et -1 (ou -1 et 1).

Tr1<-c(1,-2,1)
Tr2<-c(1,0,-1)
sex1<-c(1,-1)

Une fois que ces coefficients de contrastes ont été créés, on peut créer le vecteur du produit de ces deux vecteurs. Par exemple, on peut obtenir le produit du premier vecteur des coefficients de contrastes de la variable traitement avec le vecteur des coefficients de contrastes de la variable ‘sexe’ de la manière suivante :

c(Tr1%o%sex1)

## [1]  1 -2  1 -1  2 -1

Il ne reste plus qu’à ajouter ces produits des coefficients de contrastes à une matrice. Par souci de simplicité et pour augmenter la lisibilité, on va l’ajouter à ‘mes.coef’ pour illustrer la manière dont les contrastes sont organisés par rapport aux combinaisons de modalités.

mes.coef$cont1<-c(Tr1%o%sex1)
mes.coef$cont2<-c(Tr2%o%sex1)
mes.coef

##   Treatment    sex c.1 c.2 cont1 cont2
## 1   Control Female  -1  -1     1     1
## 2      High Female   1   0    -2     0
## 3       Low Female   0   1     1    -1
## 4   Control   Male   1   1    -1    -1
## 5      High   Male  -1   0     2     0
## 6       Low   Male   0  -1    -1     1

A présent, que la matrice de contrastes est prête, on peut utiliser la fonction contrast sur les colonnes 3 à 6 de “mes.coef”.

contrast(emmeans.out,mes.coef[,3:6] )

##  contrast estimate    SE  df t.ratio p.value
##  c.1        0.1100 0.131 309  0.841  0.4011 
##  c.2        0.0841 0.105 303  0.801  0.4236 
##  cont1     -0.1359 0.241 311 -0.563  0.5740 
##  cont2     -0.0841 0.105 303 -0.801  0.4236

5.8. Tester l’homogénéité des variances

Lorsque nous nous sommes intéressés aux statistiques descriptives, nous avons réalisé un graphique (Figure 3). Le graphique que nous avions réalisé semblait indiquer que les variances étaient plutôt similaires entre les deux groupes de traitement, mais que la variance était plus grande pour le groupe de contrôle.

Nous pouvons spécifier cette structure en précisant l’argument “weights” de la fonction lme. Par défaut, l’argument “weights” est nul, ce qui signifie que les variances sont considérées comme étant homogènes. Lorsqu’on veut préciser que les variances sont hétérogènes, il faut utiliser une formule impliquant un (ou plusieurs) facteurs fixes. Dans notre cas, l’effet fixe qui sera introduit sera l’effet du traitement.

En l’occurrence, l’hypothèse est que les variances de l’effet aléatoire de la portée ne sont pas homogènes en fonction du traitement qui a été reçu. Comme l’hypothèse est sur un facteur aléatoire, l’algorithme qui doit être utilisé est le maximum de vraisemblance restreint.

 Modele.lme.var <- lme(weight ~ Treatment + sex + Lsize, random = ~ 1 |Litter ,data = rat, method = "REML",
                     weights = varIdent(form = ~1 | Treatment))
VarCorr( Modele.lme.var)

## Litter = pdLogChol(1) 
##             Variance  StdDev   
## (Intercept) 0.0983841 0.3136624
## Residual    0.2646360 0.5144278

summary( Modele.lme.var)$modelStruct$varStruct 

## Variance function structure of class varIdent representing
##   Control       Low      High 
## 1.0000000 0.5643872 0.6354988

Dans le modèle, on constate que l’écart-type résiduel est de 0.514. Pour obtenir l’écart-type de chacun des groupes, il faut multiplier l’écart-type résiduel par les coefficients affichés dans la section “Variance function structure of class varIdent representing”. Par défaut, la modalité qui arrive en premier par ordre alphabétique a un coefficient de 1. Ainsi, pour le groupe de contrôle, l’écart-type vaut 0.514. Pour le groupe ‘Low’, l’écart-type résiduel vaut \(0.514 \times 0.565\), ce qui donne \(0.29041\). De la même manière, on peut calculer l’écart-type résiduel pour le groupe ‘High’ par le calcul suivant \(0.514 \times 0.635\), ce qui donne \(0.32639\).

Ces résultats nous incitent à penser que les variances ne sont pas homogènes entre les 3 groupes. Les modèles linéaires mixtes permettent de tester si la vraisemblance du modèle est accrue lorsqu’on considère que les variances ne sont pas homogènes.

# Modèle avec l'algorithe REML où les variances ne sont pas considérées comme hétérogènes 
 Modele.lme <- lme(weight ~ Treatment + sex + Lsize , random = ~ 1 |Litter ,data = rat, method = "REML")
anova(Modele.lme,Modele.lme.var)

##                Model df      AIC      BIC    logLik   Test  L.Ratio
## Modele.lme         1  7 414.5829 440.8952 -200.2914                
## Modele.lme.var     2  9 376.7481 410.5782 -179.3740 1 vs 2 41.83483
##                p-value
## Modele.lme            
## Modele.lme.var  <.0001

Avec un ratio de vraisemblance qui est égal à \(41.8348254\) et une valeur de probabilité qui est inférieure à 0.001, cela signifie que modéliser des variances hétérogènes améliore de manière significative le maximum de vraisemblance.

Néanmoins, on pourrait penser que les variances ne sont pas hétérogènes entre tous les groupes. Elles pourraient être homogènes pour les groupes ‘Traitements’ et se différencier du groupe ‘Contrôle’. Pour avoir une illustration sur la manière de modéliser l’hétérogénéité des variances de cette manière, vous pouvez reporter au complément à cette formation.

6. Les facteurs en mesure répétées

Pour illustrer les facteurs en mesure répétée, nous pouvons nous appuyez sur les résultats de Douglas et al. (2004). Dans cet exemple, les auteurs ont mesuré l’activité cérébrale suite à un traitement médicamenteux chez de rats au niveau de 3 régions particulière. Plus précisément, leur but était d’examiner l’activation de nucléotide dans noyaux cérébraux différents (les régions cérébrales). Ils ont mesuré cette activité après l’injection d’un solution saline et l’ont comparé à l’activation dans la même région après l’injection de carbachol. Cette activité cébébrale est produite par autoradiographie.

Le jeu de données comprend les variables suivantes :

• Animal : l’identifiant de l’animal
• Le traitement : Carbachol vs. solution saline (baal)
• La localisation cérébrale : BST, LS et VDB

Cependant, à l’heure actuelle, il est présenté en format large et il est nécessaire de le faire passer en format long.

## Classes 'tbl_df', 'tbl' and 'data.frame':    5 obs. of  7 variables:
##  $ animal   : chr  "R111097" "R111397" "R100797" "R100997" ...
##  $ Carb,BST : num  372 493 665 515 589
##  $ Carb,LS  : num  302 356 587 438 494
##  $ Carb,VDB : num  450 460 727 604 621
##  $ Basal,BST: num  366 376 458 480 463
##  $ Basal,LS : num  199 205 245 261 278
##  $ Basal,VDB: num  187 179 237 196 262

6.1. Le formatage des données : le format long et le format large.

Habituellement, on organise un jeu de données de la manière suivante : une colonne par variable et une ligne participant. Les choses se compliquent lorsque nous avons des variables en mesure répétées. Deux modes de présentation existent : le format large et le format long.

Dans le format large, pour les facteurs en mesure répétées, une colonne représente une modalité de la variable à mesure répétée (ou une combinaison de modalités s’il y a plusieurs variables). Dans l’exemple ci-dessous, il y a deux variables : séance et choc. Pour préciser à quelle séance et à quelle condition de choc appartient la mesure, il faut donner deux informations dans l’entêtes de colonnes. Ceci crée de l’ambiguïté.

choc<-read.csv("choc.csv", sep=";")
str(choc)

## 'data.frame':    30 obs. of  7 variables:
##  $ groupe        : Factor w/ 2 levels "bas","haut": 2 2 2 2 2 2 2 2 2 2 ...
##  $ sean1sans_choc: int  120 123 101 112 113 108 117 119 112 124 ...
##  $ sean2sans_choc: int  92 95 94 89 88 86 85 95 92 97 ...
##  $ sean3sans_choc: int  24 30 22 24 27 28 31 25 24 27 ...
##  $ sean1choc     : int  148 120 143 142 118 143 128 151 124 138 ...
##  $ sean2choc     : int  106 118 117 93 124 118 121 104 85 114 ...
##  $ sean3choc     : int  124 122 128 131 123 110 109 122 118 110 ...

Pour faire face à ce problème, on peut utiliser un format long. Dans le format long, une colonne représente une variable. Les mesures répétées sont identifiables par le fait que l’identifiant des participants (la variable aléatoire) va être répété.

Pour passer du format large au format long, il existe plusieurs possibilités. Une possibilité, relativement économique lorque vous avez plusieurs facteurs en mesure répétée est d’utiliser la fonction ez.reshape de easieR.

Les arguments de cette fonction sont les suivants :

• data : jeu de donnée sur lequel faire l’analyse
• varying : nom des variables qui dans le format large correspondent à des variables individuelles dans le format long
• idvar : le nom des varianle qui correspondent à l’identifiant des individus. On y inclut aussi les facteurs à groupes indépendants
• times : nom des variables mesurées en format long
• IV.names : liste avec le nom des variables indépendantes créées dans le format long
• IV.levels : liste avec les niveaux pour chaque variable indépendant créée dans le format long

De manière concrète, dans notre exemple, la fonction s’utilise de la manière suivante :

choc.long<-ez.reshape(data=choc, 
           varying =list(c('sean1sans_choc','sean2sans_choc',
          'sean3sans_choc','sean1choc','sean2choc','sean3choc')),
           v.names =c('distance'),idvar =c('groupe'),
          IV.names=list('choc','seance'),
          IV.levels=list( c('sans','avec'), c('s1','s2','s3')))

## ez.reshape(data=choc, varying =list(c('sean1sans_choc','sean2sans_choc','sean3sans_choc','sean1choc','sean2choc','sean3choc')), v.names =c('distance'),idvar =c('groupe','IDeasy'),IV.names=list('choc','seance'), IV.levels=list( c('sans','avec'), c('s1','s2','s3')))

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following object is masked from 'package:car':
## 
##     recode

## The following object is masked from 'package:nlme':
## 
##     collapse

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

head(choc.long)

##   groupe IDeasy           time distance choc seance
## 1   haut     p1 sean1sans_choc      120 sans     s1
## 2   haut     p2 sean1sans_choc      123 sans     s1
## 3   haut     p3 sean1sans_choc      101 sans     s1
## 4   haut     p4 sean1sans_choc      112 sans     s1
## 5   haut     p5 sean1sans_choc      113 sans     s1
## 6   haut     p6 sean1sans_choc      108 sans     s1

Pour les personnes qui préfèrent cliquer, vous pouvez également utiliser les boîtes de dialogue en utilisant la fonction ez.reshape() sans argument.

1. Vous commencez par choisir le jeu de données s’il y en a plusieurs en mémoire.
2. Vous indiquez le nombre de variables dépendantes. En l’occurrence, il n’y en a qu’une.
3. vous choisissez les colonnes qui correspondent aux facteurs en mesure répétées. En l’occurrence, les 6 dernières colonnes (‘sean1sans_choc’,‘sean2sans_choc’,‘sean3sans_choc’,‘sean1choc’,‘sean2choc’,‘sean3choc’)
4. Vous donnez un nom à la variable dépendante. Ici, il s’agit de la distance par rapport au centre la cible. Il faut un nom simple, sans accent, sans espace, sans caractère spéciaux. Par exemple, “distance”.
5. Ensuite, il vous est demandé combien de facteurs en mesure répétée il y a. En l’occurrence, il y a la séance et la présence du choc. Donc 2.
6. Vous devez entrer le nom de la 1e variabe en mesure répétée. Dans la structure de votre base de données, le nom de la variable à entrer en premier correspond au nom de la variable dont les modalités varient le plus lentement quand on se déplace de colonnes. En l’occurrence, quand on passe de la 2e colonne (“sean1sans_choc”) à la 3e colonne (“sean2sans_choc”), la variable séance varie (on passe de la séance 1 à la séance 2) alors que la variable choc ne varie pas (on est dans la modalité ‘sans choc’). Il faut donc entrer la variable “choc” en premier.
7. Il vous est demandé le nombre de modalités. En l’occurrence, il y en a 2.
8. Dans la console (là où on tape les lignes de commandes et où les résultats s’affichent), vous devez entrer le nom des modalités. Ici, il s’agit de “sans” et “avec”, dans l’ordre d’apparition des colonnes.
9. On répéte l’opération pour la seconde variable
10. le jeu de données est passé en format long, il vous est présenté et il vous est demandé d’appuyer sur la touche entrée quand la console est sélectionnée. Vous devez le faire après avoir vérifié si les variables ‘choc’ et ‘séance’ correspondent à l’association des modalités qu’on retrouve dans la variable ‘time’.
11. Une boîte de dialogue demande la confirmation de la bonne structure du jeu de données. Si la structure est bonne, cliquez sur “oui”.

Votre jeu de données est présent en format long et dans la mémoire de R avec comme nom le nom de votre jeu de données et une extension “.long”. Ainsi, dans notre exemple, il s’agit de “choc.long”

Exercice 13 : faites passer le jeu de données du format large au format long.

Exercice 14 : Réalisez un modèle linéaire mixte avec la fonction lme en considérant qu’il n’y a pas de facteurs en mesure répétées

6.1.Spécifier les facteurs en mesure répétées

Pour spécifier les facteurs en mesures répétées sur les facteurs aléatoires, il existe au moins deux possibilités : modéliser le croisement entre un effet aléatoire et un effet fixe ou de calculer une pente aléatoire.

6.1.1. Spécifier le croisement entre un effet aléatoire et un effet fixe

Une première manière de préciser qu’une variable est en mesure répétée est de considérer que le facteur est croisé avec l’effet aléatoire. Dans ce cas, on peut lire l’effet aléatoire comme l’effet de l’animal et l’effet du traitement à l’intérieur de l’animal. Pour spécifier ce croisement, on divisera l’effet aléatoire par l’effet fixe avec lequel il y a un croisement. Pour le présenter autrement, l’effet aléatoire prendra la forme :

\[1|effet\ aléatoire /effet\ fixe\]

Dans notre exemple, on peut donc formaliser que le traitement est en mesure répétée de la manière suivante :

 modele1a<-lme(Activite~Traitement*Region, random=~1|animal/Traitement, data= ratbrain.long, method="ML")

6.1.2. Spécifier une pente aléatoire

La seconde possibilité consiste à indiquer qu’il est nécessaire d’avoir un intercept spécifique à chaque niveau de facteur fixe mais également une pente qui lui sera propre. De manière formelle, l’écriture de l’effet aléatoire prendre la forme suivante :

\[effet\ fixe\ |\ effet\ aléatoire\]

Par exemple, si on veut modéliser un intercept pour chaque modalité du niveau de la variable Traitement, on devra exécuter la commande suivante :

 modele1b<-lme(Activite~Traitement*Region, random=~Traitement|animal, data= ratbrain.long, method="ML")

6.1.3. Croisement entre l’effet fixe et l’effet aléatoire ou pente aléatoire

La différence majeure entre le modèle 1a et le modèle 1b est le nombre de paramètres estimés : dans le modèle 1b, comme on calcule une constante pour chaque niveau d’une variable indépendante, on va utilisé des degrés de liberté supplémentaires (correspondant au nombre de modalités de la variable moins -1). En revanche, dans le cas où on précise que le facteur est croisé, on n’utilise qu’un de degré de liberté supplémentaire. On peut constater ce phénomène en réalisant l’anova entre les deux modèles, et en le comparant avec la situation où il n’y a ni croisemen entre l’effet fixe et l’effet aléatoire, ni pente aléatoire :

anova(modele0, modele1a, modele1b)

##          Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## modele0      1  8 341.3376 352.5472 -162.6688                         
## modele1a     2  9 320.1284 332.7391 -151.0642 1 vs 2 23.209237  <.0001
## modele1b     3 10 312.7230 326.7349 -146.3615 2 vs 3  9.405403  0.0022

Ainsi, le modèle peut être moins parcimonieux quand on estime une pente par niveau d’une variable aléatoire. Il faut donc éviter d’ajouter des constantes par facteur si cela ne s’avère pas nécessaire. Il faut donc s’assurer de la pertinence de le faire. Pour vous aider dans cette démarche, une inspection visuelle préalable des données peut être utile. Une manère classique de représenter l’évolution des données individuelles est de faire une représentation d’un graphique en treillis avec comme particularité de faire la représentationniveau de l’effet aléatoire. On peut obtenir ce type de graphique par la fonction interaction.plot. En l’occurrence, nous allons réaliser un graphique d’interaction entre notre variable aléatoire et un facteur en mesure répétée.

 interaction.plot(ratbrain.long$Traitement, ratbrain.long$animal, ratbrain.long$Activite, las=1,
                  trace.label="identifiant du rat", xlab="Traitement", ylab="Activité cérébrale")

Dans ce graphique, on se rend compte que les différents rats n’ont pas la même ligne de base et surtout que l’impact du carbachol n’est pas le même pour les différents rats. Dans ce cas, il est utile de modéliser une pente aléatoire pour chaque rat.

Cela est confirmé par l’anova que nous avons réalisé précédemment qui montrait que la modélisation d’une pente aléatoire améliorait significativement la vraisemblance du modèle.

En termes d’interprétation, cela signifie que l’effet de la variable indépendante n’est pas le même pour tous les niveaux de l’effet aléatoire. Dans notre exemple, cela signifie que l’effet du traitement ne sera pas le même pour l’ensemble des rats.

A titre pédagogique, il est intéressant d’identifier une erreur qui aurait pu être faite : considérer l’effet de la région cérébrale comme un facteur en mesure répétée. En effet, pour les personnes habituées avec les analyses de variance, il serait assez naturel de le considérer comme tel. En l’occurrence, même si les mesures sont faites chez le même rat, ce n’est pas une mesure répétée à proprement parlé puisque la mesure répétée implique de refaire la mesure sur le même échantillon. En l’occurrence, les trois régions cérébrales sont 3 échantillons différents, dont on peut penser que les mesures soient plus proches entre elles quand elles appartiennent à un cerveau plutôt qu’à des cerveaux différents.

Une autre manière de voir les choses, c’est de réaliser une inspection visuelle des données concernant l’effet de l’activité cérébrales afin de voir si cela aurait un intérêt de modéliser ce facteur comme un facteur en mesure répétée. En l’occurrence, on observe que l’effet des régions cérébrales est pratiquement le même pour chaque rat, rendant plutôt inutile de modéliser cet effet dans notre modèle.

  interaction.plot(ratbrain.long$Region, ratbrain.long$animal, ratbrain.long$Activite, las=1,
                   trace.label="identifiant du rat", xlab="Région cérébrale", ylab="Activité cérébrale")

En effet, lorsque le modèle inclut comme facteur en mesure répété les régions, la vraisemblance n’est pas améliorée.

  modele2a<-lme(Activite~Traitement*Region, random=~1|animal/Region, data= ratbrain.long, method="ML")
  modele2b<-lme(Activite~Traitement*Region, random=~Region|animal, data= ratbrain.long, method="ML")
  anova(modele0, modele2a, modele2b)

##          Model df      AIC      BIC    logLik   Test     L.Ratio p-value
## modele0      1  8 341.3376 352.5472 -162.6688                           
## modele2a     2  9 343.3376 355.9484 -162.6688 1 vs 2 2.87045e-08  0.9999
## modele2b     3 13 351.3376 369.5532 -162.6688 2 vs 3 7.70558e-09  1.0000

6.1.4. La variable est-elle ou non en mesure répétée ?

A présent que vous savez modéliser des facteurs en mesure répétée, il faut se poser la question de savoir quels facteurs en mesure répétée il faut modéliser. Dans une anova classique, on aurait tendance à modéliser à la fois l’effet de la variable région cérébrale et l’effet du médicament comme facteur en mesure répétée. Dans le cas du modèle linéaire mixte, les choses sont un peu différentes : il faut se rappeler qu’on doit pouvoir calculer une variance sur chaque modalité de l’effet aléatoire. Si vous n’avez qu’une mesure pour chaque modalité de l’effet aléatoire, vous ne pourrez pas calculer de variance.

Cela étant dit, un autre élément doit entrer en considération : en réalité dans l’anova, les modalités des facteurs en mesure répétée doivent être predictifs les unes des autres. Connaître le score au temps 1 est prédictif du score au temps 2 et au temps 3. Cette règle est cependant souvent ignorée dans le cas de l’anova (raison pour laquelle on utilise des correction des degrés de liberté pour y faire face) Cette règle vaut aussi pour les modèles linéaires mixtes. A la différence que nous n’avons pas besoin de l’ignorer. En effet, on peut se poser la question de savoir si une modalité est prédictive d’une autre modalité. Si la réponse est non, on ne doit pas la modéliser. On peut difficilement imaginer que l’activité cérébrale à un moment T soit prédictible de l’activité des autres régions à ce moment T. Dès lors, il semble peu pertinent de le modéliser.

Pour vous en convaincre, nous allons tout de même le modéliser et comparer avec une situation où seul l’effet du traitement est pris en considération. Pour modéliser deux facteurs en mesure répétées, on comment par le facteur du plus haut niveau qu’on divise par un facteur de plus bas niveau (le traitement est inclus dans toutes les modalités de la variable Région).

  modele3a<-lme(Activite~Traitement*Region, random=~1|animal/Region/Traitement, data= ratbrain.long, method="ML")

Quand on compare ce nouveau modèle avec les précédent, on constate que le fait d’introduire la région dans le modèle détériore considérablement la vraisemblance du modèle pour revenir au niveau initial.

  anova(modele0, modele1b, modele3a)

##          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## modele0      1  8 341.3376 352.5472 -162.6688                        
## modele1b     2 10 312.7230 326.7349 -146.3615 1 vs 2 32.61464  <.0001
## modele3a     3 10 345.3376 359.3496 -162.6688

Pour illustrer de manière encore plus forte le côté non pertinent de modéliser les deux, la fonction lmer du package lme4 va renvoyer un message d’erreur si vous essayez de modéliser un effet aléatoire pour lequel il n’y a qu’une modalité. En l’occurrence, il n’y a qu’une mesure par animal au croisement entre les régions et les traitement, empêchant dès lors toute variance d’être calculée.

  modelelmer<-try(lmer(Activite~Traitement*Region+(1|animal/Region/Traitement), data= ratbrain.long, REML=F))

## Error : number of levels of each grouping factor must be < number of observations

modelelmer

## [1] "Error : number of levels of each grouping factor must be < number of observations\n"
## attr(,"class")
## [1] "try-error"
## attr(,"condition")
## <simpleError: number of levels of each grouping factor must be < number of observations>

Exercice 15 : Modéliser des variances hétérogènes sur la variable “Traitement” et tester si cela améliore le modèle 1b

La procédure pour préciser les facteurs en mesure répétées est la même pour la fonction lmer que pour la fonction lme, vous pouvez essayer la modélisation avec la fonction lmer.

Exercice 16 : Obtenez les contrastes d’interaction sur le modèle 1b et faites-en une représentation graphique

6.2. Les modèles avec pentes aléatoires pour les données longitudinales ou analyse de courbe de croissance

6.2.1. Enoncé de travail

Dans cet exemple, des auteurs s’intéressent au développement social d’enfants autistes. Pour chaque enfant, un score de socialisation a été mesuré à 2, 3, 5, 9 et 13 ans. On va prendre compte un modèle qui permet à chaque enfant d’avoir un coefficient qui varie de manière aléatoire permettant ainsi de décrire un trajectoire individuelle.

Les modèles avec des coefficients aléatoires sont souvent utilisés dans le cas d’études longitudinales lorsque le chercheur s’intéresse dans la modélisation des effets du temps, ainsi que de la modélisation de covariables qui varie avec le temps.

Dans le cas particulier de la psychologie du développement, on fait souvent référence aux modèles avec coefficients aléatoires comme des modèles des courbes de croissance.

Pour illustrer cette thématique, nous allons utiliser des données issues de l’université de Michigan (Anderson et al., 2009). Cette étude a été réalisée de manière longitudinale sur 214 enfants. Ces enfants ont été divisés en trois catégories diagnostiques lorsqu’ils avaient 2 ans : autisme, trouble envahissant du développement (TED -PDD en anglais) et enfant au développement typique. Le jeu de données à votre disposition comprendre 158 enfants du spectre autistique (y compris les TED). L’étude a été pensée pour avoir des informations à chacun des âges (2, 3, 5, 9 et 13 ans), bien que tous les enfants n’aient pas été vus à tous les âges. Un des objectifs de l’étude était d’évaluer l’influence relative de la catégorie diagnostique initiale, du niveau langagier à 2 ans et d’autres covariables sur les trajectoires développementales de la socialisation de ces enfants.

Le développement social était évalué par le Vineland Adaptive Behavior Interview (VSAE). Ce questionnaire évalue les relations interpersonnelles, le temps de jeu, les capacités de coping. Le niveau langagier initial a été évalué par le Sequenced Inventory of Communication Developement (SICD). Les enfants étaient placés dans un des trois groupes (sicdegp) sur la base de leur score initial au SICD à 2 ans.

autisme<-read_excel("autism.xlsx", sheet="autism")
str(autisme)

## Classes 'tbl_df', 'tbl' and 'data.frame':    612 obs. of  4 variables:
##  $ age    : num  2 3 5 9 13 2 3 5 9 13 ...
##  $ vsae   : num  6 7 18 25 27 17 18 12 18 24 ...
##  $ sicdegp: num  3 3 3 3 3 3 3 3 3 3 ...
##  $ childid: num  1 1 1 1 1 3 3 3 3 3 ...

Dans ce jeu de données, les variables correspondent à : - age : l’âge de l’enfant - vsae : la variable dépendante, le score de développement social - scicdegp : le niveau langagier (1 = faible, 2 = moyen, 3 = haut)

Exercice 17 : dans une colonne appelée “lang”, recoder les niveaux langagier en faible, moyen et haut et recoder l’âge en variable catégorielle dans une autre colonne.

Comme les enfants ont été testés la 1e fois à 2 ans, la valeur prédite VSAE doit être celle qu’on obtient à 2 ans. Cependant, la constante, par défaut considère que c’est 0. Pour simplifier donc l’interprétation, il suffit de retirer deux ans à la variable âge de sorte à ce que la constante à l’origine correspondent au VSAE prédit à deux ans. On réaliser cette opération de la manière suivante :

autisme$age.2<-autisme$age-2

6.2.2. Le graphique en treillis

Nous avons vu précédemment qu’on pouvait obtenir un graphique en treillis avec la fonction interaction.plot afin de voir l’évolution temporelle spécifique à chaque niveau (ici enfant autiste) de l’effet aléatoire. Dans cette section, nous allons aborder d’autres manières d’obtenir ce graphique de sorte à gagner en flexibilité. Les deux packaes qui seront utilisés sont ggplot2 et lattice. Une description sera faite pour chaque package, à vous de choisir la manière qui vous convient le mieux pour vos graphiques.

Exercice 18 : Faites une représentation graphique (graphique en treillis) qui permet de voir l’évolution individuelle de chaque enfant en fonction de l’âge comme dans la Figure ci-dessous.

Sur la base de ce graphique, on identifie : 1) qu’il y a eu une perte expérimentale : on a les mesures à 2 et à 3 ans pour certains enfants mais pas après 4 ans. Pour d’autres les mesures commencent à 5 ou à 9 ans. On peut objectiver ce constater en regardant le nombre d’observations pour chacun des âges grâce à la fonction summary :

summary(autisme$age.f)

##   2   3   5   9  13 
## 156 150  91 120  95

2. que tous les enfants ne vont pas progresser de manière similaire ce qui justifie le fait de modéliser une pente aléatoire.

A ces deux premiers constats, une analyse du plan expérimental permet de voir que le plan n’est pas équilibré sur les différents niveaux de langage grâce à la fonction table.

table(autisme$age.f, autisme$lang)

##     
##      faible moyen eleve
##   2      50    66    40
##   3      48    64    38
##   5      29    36    26
##   9      37    48    35
##   13     28    41    26

Si le graphique présenté ci-dessus répond aux attentes minimales, il manque de finesse : on ne sait par exemple pas si ces évolutions différentielles s’observent pour chacun des groupes d’enfants. Ainsi, il serait utile de faire une représentation d’un graphique en treillis avec comme particularité de faire la représentation par groupe.

Il existe plusieurs manière de faire cette représentation. La première est d’utiliser le package ggplot2. De manière succincte, le graphique s’obtient en utilisant une grammaire. La fonction ggplot cherche un jeu de données (autisme) et les variables qui vont servir pour l’esthétisme du graphique (aes). En d’autre termes, aes est une fonction. Dans notre exemple, il y a 4 arguments qui vont être utiliser : x est la variable sur l’abscisse, y sur l’ordonnée, coulour est utilisé pour utiliser des couleurs différentes pour chaque modalité de la variable indiquée et group permet de respécifier à quel groupe appartient chaque observation. En l’occurrence, il faut répéter group et colour pour que le graphique affiche l’évolution en fonction des différents groupes d’âge. Dans le cas contraire, par défaut, il fera un ligne verticale par groupe d’âge. Ensuite, il faut préciser le type de graphique (la géométrie) désiré, c’est le “geom”. Pour un graphique en ligne, cela s’écrit par la fonction geom_line(). On termine en précisant qu’on voudrait avoir cette représentation graphique en ayant un graphique par groupe langagier. Cela est possible grâce à la fonction facet_grid. Cette fonction s’utilise avec un tilde. Les variables à gauche du tilde crée des graphiques sur un axe vertical et les variables à droite sur un axe horizontale. Pour mettre plusieurs variables sur un même axe, on utilise le symbole “+”.

p<-ggplot(autisme, aes(x=age.f, y=vsae, colour=childid, group=childid))
p<-p+geom_line()
p<-p+facet_grid(~lang)
p

Une alternative est d’utiliser la fonction plot.trellis du package lattice. Pour l’utiliser, on commence par créer un objet “groupedData”. On précise par une formule la variable dépendante, la variable indépendante d’intérêt et le facteur aléatoire. Pour avoir un graphique pour chaque niveau de langage, on utilise l’argument “outer” sur un principe similaire à ce qu’on observe pour ggplot2.

autism.gl <- groupedData(vsae ~ age | childid,
outer = ~ lang, data = autisme)
plot(autism.gl, display = "childid", outer = TRUE,
key = F, xlab = "Age", ylab = "VSAE",
main = "Données individuelle en fonction du niveau langagier ")

On peut également faire une représentation graphique pour laquelle nous représentatons la moyenne à chaque groupe d’âge, ce qui permet de voir l’évolution des compétences sociales pour chacun des groupes. La construction du graphique s’appuie sur le même principe que précédemment : puisqu’on veut la moyenne du niveau langagier le groupedData doit inclure la variable “lang” comme facteur après la barre verticale et qu’on affiche la moyenne du niveau langagier en annonçant que ce qui doit être affiché (display) est “lang” et pas “childid”

Dans les trois graphiques, l’information importante à observer est que, s’il y a une variabilité développementale sur le score du VSAE, cette variabilité est moins marquée pour le groupe avec un faible niveau langagier que pour les deux autres groupes.

6.2.3. Modélisation de données longitudinales.

6.2.3.1. L’intérêt de la fonction groupedData

Pour commencer, on va réaliser un jeu de données “regroupé” en utilisant la fonction groupedData. Cela va permettre de spécifier la nature hiérarchique du jeu de données. Les arguments de cette fonction indique que VSAE est la variable dépendante, age.2 est la première covariable et que “childid” définit les “groupes” d’observations auxquels sont associés les effets aléatoires.

autisme.gr<-groupedData(vsae~age.2|childid, data= autisme)
head(autisme.gr)

## Grouped Data: vsae ~ age.2 | childid
##   age vsae sicdegp childid  lang age.f age.2
## 1   2    6       3       1 eleve     2     0
## 2   3    7       3       1 eleve     3     1
## 3   5   18       3       1 eleve     5     3
## 4   9   25       3       1 eleve     9     7
## 5  13   27       3       1 eleve    13    11
## 6   2   17       3       3 eleve     2     0

Pour le formuler autrement, nous venns d’associer un modèle à notre jeu de données. Cela facilite des opérations ultérieures qu’on voudrait réaliser. Par exemple, nous pouvons à présent réaliser un modèle linéaire mixte ainsi :

 modele.test.groupeddata<-lme(autisme.gr, na.action=na.omit)
summary(modele.test.groupeddata)

## Linear mixed-effects model fit by REML
##  Data: autisme.gr 
##        AIC     BIC  logLik
##   4767.199 4793.66 -2377.6
## 
## Random effects:
##  Formula: ~age.2 | childid
##  Structure: General positive-definite
##             StdDev    Corr  
## (Intercept) 0.3899919 (Intr)
## age.2       4.3547658 1     
## Residual    7.9132485       
## 
## Fixed effects: vsae ~ age.2 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 9.517918 0.4571426 451 20.82046       0
## age.2       4.416744 0.3759710 451 11.74757       0
##  Correlation: 
##       (Intr)
## age.2 -0.144
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.74904702 -0.37191669 -0.06259569  0.31866362  5.45892701 
## 
## Number of Observations: 610
## Number of Groups: 158

On peut également réaliser des graphiques assez aisément, directement sur le jeu de données :

 plot(autisme.gr)

En résumé, la fonction groupedData n’est pas obligatoire, mais peut faciliter certains traitement ultérieurs.

6.2.3.1.Spécifier des relations non linéaires

Nous allons à présent formaliser notre modèle en incluant :

• une relation linéaire entre l’âge et le VSAE
• une relation quadratique entre l’âge et le VSAE
• une relation linéraire entre le niveau langagier et le VSAE
• une interaction entre âge et niveau langagier
• une interaction entre la relation quadratique de l’âge et le niveau langagier

Pour spécifier une relation quadratique qui soit considérée comme telle, il faut expliquer à R qu’il doit intérpréter le carré. Pour cela, on va utiliser la fonction I, qui signifie “asIs”.

# remarquer que comme on utilise le jeu de données "autisme.gr", des informations 
# peuvent être omises dans le modèle le rendant plus lisible ici. 
try(autisme1<-lme(vsae~age.2+I(age.2^2)+lang+
                  age.2:lang+I(age.2^2):lang,
              random=~age.2+I(age.2^2)-1, method="REML", data= autisme.gr, na.action=na.omit))

Comme le modèle ne peut pas converger avec le calcul d’un constante pour chaque enfant, on va retirer ce paramètre des effets aléatoires.

autisme2<-lme(vsae~age.2+I(age.2^2)+lang+
                  age.2:lang+I(age.2^2):lang,
              random=~age.2+I(age.2^2)-1, method="REML", 
              data= autisme.gr, na.action=na.omit)

On peut se demander si la relation quadratique pour l’effet aléatoire est utile ici.

Exercice 19 : Testez si l’effet spécifique de la relation quadratique est utile.Voici les résultats que vous devriez obtenir :

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme2      1 13 4641.276 4698.457 -2307.638                        
## autisme2a     2 11 4721.203 4769.587 -2349.601 1 vs 2 83.92683  <.0001

Nous allons à présenter essayer de trouver le modèle le plus parcimonieux pour tester l’effet de la progression sur le score de VSAE. Pour cela, nous allons voir si l’effet d’interaction entre l’âge au carré (la relation quadratique) et le niveau langagier est utile. Remarquez que, pour m’assurer du modèle qui est testé, la fonction formula est utilisée afin de connaître les termes du modèle.

## vsae ~ age.2 + I(age.2^2) + lang + age.2:lang + I(age.2^2):lang
## <environment: 0x0000000031969eb0>

## vsae ~ age.2 + I(age.2^2) + lang + age.2:lang
## <environment: 0x000000003e2031d8>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3      1 13 4636.444 4693.819 -2305.222                        
## autisme3a     2 11 4634.314 4682.862 -2306.157 1 vs 2 1.869704  0.3926

En l’occurrence, ce n’est pas le cas, le modèle 3a est ici plus parcimonieux.

Nous allons à présenter tester si le second terme d’interaction (entre le langage et l’âge) a un apport significatif.

## vsae ~ age.2 + I(age.2^2) + lang
## <environment: 0x000000003e2031d8>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3b     2  9 4653.696 4693.417 -2317.848 1 vs 2 23.38232  <.0001

En l’occurrence, l’interaction est significative, il est utile de la maintenir dans le modèle.

Nous allons à présenter tester si 1) l’effet du niveau langagier a un apport significatif au modèle.

## vsae ~ age.2 + I(age.2^2) + age.2:lang
## <environment: 0x000000003e2031d8>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3c     2  9 4653.929 4693.650 -2317.964 1 vs 2 23.61505  <.0001

2. l’effet de la relation quadratique de l’âge a un effet significatif

## vsae ~ age.2 + lang + age.2:lang
## <environment: 0x000000003e2031d8>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3d     2 10 4638.766 4682.900 -2309.383 1 vs 2 6.451775  0.0111

3. l’effet de l’âge a un effet significatif

## vsae ~ I(age.2^2) + lang + lang:age.2
## <environment: 0x000000003e2031d8>

##           Model df      AIC      BIC    logLik
## autisme3a     1 11 4634.314 4682.862 -2306.157
## autisme3e     2 11 4634.314 4682.862 -2306.157

Autant, le niveau langagier et la relation quadratique à l’âge ont un effet significatif, autant il n’y a pas de différence entre le modèle avec l’âge ou sans l’âge, à la condition d’avoir la relation quadratique. Cela s’explique par le fait que, pour la relation quadratique, l’effet linéaire est implicitement testée ici. Néanmoins la table de l’anova révèle que tous les effets inclus dans le modèle 3a sont significatifs.

##             numDF denDF  F-value p-value
## (Intercept)     1   448 947.8164  <.0001
## age.2           1   448 146.4845  <.0001
## I(age.2^2)      1   448   6.2236   0.013
## lang            2   155  23.0849  <.0001
## age.2:lang      2   448  12.6714  <.0001

6.2.3.4.Spécifier la structure de la matrice devariance et covariance

Une des particularité des modèles linéaires mixtes est qu’on peut spécifier la structure de la matrice de corrélation intragroupe. Pour le formuler autrement, on peut imaginer qu’un enfant à 12 ans ressemble plus à ce qu’il était lorsqu’il avait 10 ans que lorsqu’il en avait 4.

De manière théorique, il est possible de créer notre propre fonction spécifiant l’organisation de la matrice de corrélation, mais il existe une multitude de structure de matrice de corrélations qui sont proposées dans le package nlme. C’est ce que Bates et Pinheiro appellent les corStruct, qui ont deux arguments principaux : “value” et “form”. Le premier argument spécifie la valeur des paramètres de corrélations et le second est une formule à un côté qui prend la forme d’une variable temporelle (ou plus largement positionnelle) et l’effet aléatoire.

\[form = ~ age\ |\ enfant\]

Avant de pouvoir modéliser la matrice de corrélation intragroupe, il est nécessaire de connaître les différentes structures de matrices de corrélations afin de pouvoir choisir la structure la plus adaptée.

• corCompSymm : cette fonction modélise une matrice de corrélation respectant la sphéricité. Il s’agit de la structure la plus simple. Elle suppose que toutes les corrélations sur les erreurs intra-groupe soient égales. Cela signifie que la corrélation entre deux moments ne dépend pas de leur position relative. En d’autres termes, la corrélation entre 2 ans et 10 ans et la même qu’entre 6 ans et 10 ans. Par exemple :

\[corCompSymm( value = 0.3, form = ~ 1 | childid )\] signifie qu’il y a sphéricité de la matrice de covariance avec un coefficient de corrélation intraclasse de 0.3.

De manière concrète, une matrice de corrélation sphérique est une matrice dont la diagonale vaut 1 et les valeurs en dehors de la diagonale sont toutes identiques.

cs1 <- corCompSymm()
corMatrix(cs1, covariate = 1:3)

##      [,1] [,2] [,3]
## [1,]  1.0  0.5  0.5
## [2,]  0.5  1.0  0.5
## [3,]  0.5  0.5  1.0

Cette structure de corrélation est souvent trop simple pour des applications pratiques impliquant des séries temporelles. En général, il est plus réaliste de considérer un modèle dans lequel la corrélation entre deux observations diminue, en valeur absolue, avec leur distance.

En revanche, cette structure est utile lorsque toutes les observations sont collectées au même moment. Dans le cas contraire, il est préférable d’utiliser une structure de corrélations autorégressive et moyenne-mobile.

• Générale : il s’agit de l’extrême opposé à la sphéricité. Chaque corrélation est représentées par un paramètre différent. Comme le nombre de paramètres augmente de manière quadratique avec le nombre maximum d’observations à l’intérieur d’un groupe, cette structure aboutit souvent un sur-paramétrisation du modèle. Elle est donc rarement adaptée.

• corAR1 : il s’agit d’une structure de corrélation autorégressive d’ordre 1. Il s’agit du modèle autorégressif le plus simple et l’un des plus utiles. Les valeurs des corrélations diminuent en valeur absolue de manière exponentielle avec le décalage.

cs1 <- corAR1(0.8)
corMatrix(cs1, covariate = 1:5)

##        [,1]  [,2] [,3]  [,4]   [,5]
## [1,] 1.0000 0.800 0.64 0.512 0.4096
## [2,] 0.8000 1.000 0.80 0.640 0.5120
## [3,] 0.6400 0.800 1.00 0.800 0.6400
## [4,] 0.5120 0.640 0.80 1.000 0.8000
## [5,] 0.4096 0.512 0.64 0.800 1.0000

Dans la matrice ci-dessus, la corrélation entre deux moments successifs est de 0.80. Comme elle est d’ordre 1, s’il y a un décalage de 1 entre deux moments, la corrélation vaut \[0.8\times0.8=0.64\] Si le décalage est de 2, la corrélation vaut :

\[0.64\times0.8=0.512\] et ainsi de suite.

Lorsque cette structure de corrélations est généralisée à un temps continu, on parle de structure autorégressive continue, CAR(1).

La structure de corrélation autorégressive peut être d’un ordre supérieur à 1. L’ordre dans la structure de la matrice de corrélation représente le nombre d’observations passées inclues dans l’estimation de la corrélation.

• ARMA : il s’agit d’une structure de corrélation autorégressive avec moyenne mobile. Elle est obtenue en combinant un modèle autorégressif and un modèle avec moyenne mobile. Il y a donc deux paramètres sur ce type de structure : le premier correspond à l’ordre du modèle autorégressif (p) et le second au nombre de paramètres pour la moyenne mobile (q).

Ainsi, si vous avez compris le principe de AR1, si p vaut 1 et q vaut 0, ce modèle revient à faire AR1 ; si p vaut 0 et q>0, cela revient à faire un modèle basé sur la moyenne mobile uniquement. Ci-dessous, nous allons illustrer une structure de corrélation ARMA avec comme paramètre p =1 et q = 1.

cs2ARMA <- corARMA(c(0.8, 0.4), form = ~ 1 | childid, p=1, q=1)
cs2ARMA <- Initialize(cs2ARMA, data = autisme)
corMatrix(cs2ARMA)[[1]]

##         [,1]   [,2]  [,3]   [,4]    [,5]
## [1,] 1.00000 0.8800 0.704 0.5632 0.45056
## [2,] 0.88000 1.0000 0.880 0.7040 0.56320
## [3,] 0.70400 0.8800 1.000 0.8800 0.70400
## [4,] 0.56320 0.7040 0.880 1.0000 0.88000
## [5,] 0.45056 0.5632 0.704 0.8800 1.00000

Dans cette structure, la corrélation entre le T1 et le T2 vaut 0.88. Lorsqu’on séloigne du T1 et qu’on fait la corrélation entre le T1 et le T3, la valeur obtenue à la fois au travers du lag entre le T1 et le T2 et de la diminution de la constante de la moyenne mobile :

# entre T1 et 3
t1t3<-0.88*0.88 - (0.88)^2*0.09090909
t1t3

## [1] 0.704

# entre T1 et T4 
t1t4<-t1t3*0.88-t1t3*0.88*0.09090909
t1t4

## [1] 0.5632

# entre T1 et T5 
t1t5<-t1t4*0.88-t1t4*0.88*0.09090909
t1t5

## [1] 0.45056

• Les structures de corrélations spatiales : ces structures pourraient prendre du sens sur les données d’imagerie. Néanmoins, ces classes de structure dépassent le cadre de ce chapitre. Pour en savoir plus sur ce type de structure, on peut se référer au livre de Bates et Pinheiros (2000)

Pour connaître la structure de la matrice de corrélation qui convient le mieux à notre jeu de données, on peut utiliser la fonction d’autocorrélation ACF.

ACF( autisme3a )

##   lag         ACF
## 1   0  1.00000000
## 2   1 -0.39796730
## 3   2  0.08670546
## 4   3  0.06161871
## 5   4  0.02336154

plot(ACF(autisme3a, maxLag = 10), alpha = 0.01)

Les valeurs de la fonction d’autocorrélation indique que les valeurs des corrélations (en valeur absolue) diminuent de manière exponentielle à chaque fois que l’écart entre les délais augmente. Cette structure suggère donc qu’une structure du type autorégressif d’ordre 1 devrait permettre de rendre compte de la structure des données.

En l’occurrence, nous allons utiliser la fonction corAR1 pour modéliser l’évolution des enfants au travers du temps.

Pour spécifier la manière dont un enfant évolue, on peut utiliser l’argument correlation.

autisme4<-update(autisme3a, correlation=corAR1(form=~1|childid))

Pour savoir si modéliser la matrice de corrélation améliore la vraisemblance du modèle, on peut utiliser la fonction anova

anova(autisme3a, autisme4)

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme4      2 12 4622.090 4675.051 -2299.045 1 vs 2 14.22395   2e-04

En l’occurrence, le modèle est significativement amélioré quand on ajoute un terme spécifiant l’organisation de la matrice de variance/covariance.

A présent, lorsqu’on demande le résumé du modèle, une nouvelle information apparaît le “Phi”. Il s’agit du coefficient de corrélation entre les données à un moment T et les données au moment T+1.

summary(autisme4)

## Linear mixed-effects model fit by maximum likelihood
##  Data: autisme.gr 
##       AIC      BIC    logLik
##   4622.09 4675.051 -2299.045
## 
## Random effects:
##  Formula: ~age.2 + I(age.2^2) - 1 | childid
##  Structure: General positive-definite, Log-Cholesky parametrization
##            StdDev    Corr  
## age.2      3.9190198 age.2 
## I(age.2^2) 0.3467991 -0.335
## Residual   5.9614719       
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | childid 
##  Parameter estimate(s):
##       Phi 
## -0.238343 
## Fixed effects: vsae ~ age.2 + I(age.2^2) + lang + age.2:lang 
##                    Value Std.Error  DF   t-value p-value
## (Intercept)     8.497279 0.6147969 448 13.821277  0.0000
## age.2           2.188839 0.6485547 448  3.374950  0.0008
## I(age.2^2)      0.094315 0.0419958 448  2.245826  0.0252
## langmoyen       1.337420 0.7938480 155  1.684731  0.0941
## langeleve       4.870303 0.8928414 155  5.454835  0.0000
## age.2:langmoyen 0.547753 0.7948247 448  0.689150  0.4911
## age.2:langeleve 4.142840 0.8785408 448  4.715592  0.0000
##  Correlation: 
##                 (Intr) age.2  I(.2^2 lngmyn langlv ag.2:lngm
## age.2           -0.363                                      
## I(age.2^2)       0.217 -0.382                               
## langmoyen       -0.738  0.217  0.001                        
## langeleve       -0.655  0.191  0.004  0.508                 
## age.2:langmoyen  0.227 -0.694 -0.007 -0.313 -0.157          
## age.2:langeleve  0.207 -0.630  0.000 -0.160 -0.300  0.514   
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -4.88036257 -0.41807179 -0.05113287  0.33272466  6.77413434 
## 
## Number of Observations: 610
## Number of Groups: 158

On peut obtenir l’intervalle de confiance sur l’estimation de cette corrélaton à l’aide de la fonction intervals

intervals(autisme4)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                       lower       est.     upper
## (Intercept)      7.29598739 8.49727883 9.6985703
## age.2            0.92158642 2.18883926 3.4560921
## I(age.2^2)       0.01225695 0.09431518 0.1763734
## langmoyen       -0.22171334 1.33742029 2.8965539
## langeleve        3.11674391 4.87030258 6.6238612
## age.2:langmoyen -1.00530630 0.54775311 2.1008125
## age.2:langeleve  2.42620227 4.14283999 5.8594777
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: childid 
##                            lower       est.      upper
## sd(age.2)              3.3214136  3.9190198  4.6241505
## sd(I(age.2^2))         0.2832283  0.3467991  0.4246385
## cor(age.2,I(age.2^2)) -0.5307028 -0.3346167 -0.1045298
## 
##  Correlation structure:
##          lower      est.      upper
## Phi -0.3537541 -0.238343 -0.1157833
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##    lower     est.    upper 
## 5.520657 5.961472 6.437485

7. Plusieurs effets aléatories

8. Les modèles linéaires mixtes généralisés

Synthèse des packages et fonctions utilisées

Package	Nom de la fonction	Description
base	anova	Donne la table de l’anova
base	c	Créer un vecteur
base	cbind	Combiner des colonnes
base	coef	Donne les coefficients d’un modèle
base	contr.helmert.	Matrice de contrastes
base	contr.poly	Matrice de contrastes
base	contr.treatment	Matrice de contrastes
base	contrasts	Permet d’associer des contrastes à un facteur
base	data.frame	Créer un jeu de données
base	factor	Transformer une variable en variable catégorielle
base	getwd	Obtenir le répertoire de travail
base	head	Obtenir les premières lignes d’un jeu de données
base	install.packages	Installer un package
base	library	Charger un package
base	lm	Modéliser un modèle linéaire classique
base	matrix	Créer une matrice
base	mean	Moyenne
base	nrow	Obtenir le nombre de lignes d’un jeu de données
base	plot	Fait un graphique (la sortie dépend des données)
base	rbind	Combiner des lignes
base	require	Charger un package
base	resid	Donne les résidus bruts
base	residuals	Permet d’obtenir les résidus, notamment standardisés
base	sd	Ecart-type
base	setwd	Choisir le répertoire de travail
base	shapiro.test	Test de Shapiro-Wilk
base	summary	Donne la table des t et des coefficients de régression et MLM
base	var	Obtenir la matrice de variance et covariance
base	which.max	Obtenir la ligne contenant la valeur maximale
emmeans	contrast	Fait des contrastes élaborés
emmeans	emm_options	Permet de spécifier les options dans emmeans
emmeans	emmeans	Crée un objet compatible pour les contrastes dans emmeans
emmeans	get_emm_option	Permet de connaître les options dans emmeans
emmeans	pairs	Fait les comparaisons 2 à 2 dans emmeans
emmeans	plot	Fait une représentation graphique d’un objet emmeans
ggplot2	ggplot	Réaliser des graphiques
HLMdiag	case_delete	Déterminer l’impact de la suppression des observations
HLMdiag	diagnostics	Calculer les indicateurs statistiques d’influence
HLMdiag	dotplot_diag	Permet d’avoir une représentation graphique des indicateurs statistiques
lattice	dotplot	Permet de représenter (notamment) le facteur aléatoire
lattice	splom	Permet de déterminer l’influence des paramètres les uns par rapport aux autres.
lattice	xyplot	Permet de déterminer si les variances suivent une distribution normale
lme4	lmer	Permet de tester un modèle linéaire mixte
lme4	profile	Permet de déterminer l’intervalle de confiance sur les paramètres estimés (en particulier l’effet aléatoire)
lme4	ranef	Permet d’obtenir le BLUP
lmerTest	lmerTest::anova	Permet d’obtenir la table de l’anova (type I et III) pour lme4
lmerTest	ranova	Permet de tester si l’effet aléatoire améliore la vraisemblance
lmerTest	step	Donne les facteurs qui ont un apport significatif
MuMIn	r.squaredLR	Fournit un pseudoR² de Nagelkerke
MuMIn	squaredGLMM	Fournit le coefficient de détermination (R²)
nlme	gls	Permet de tester un modèle basé sur les moindres carrés généralisés
nlme	lme	Permet de tester un modèle linéaire mixte
nlme	lmeControl	Permet de paramétrer les valeurs pour estimer le MV
nlme	random.effects	Permet d’obtenir le BLUP
nlme/lme4	VarCorr	Obtenir les informations sur les variances de l’effet aléatoire et des résidus
outliers	grubbs.test	Identifier la présence de valeurs influentes
phia	interactionMeans	Permet d’avoir les moyennes et erreurs-types ajustées (intéressant quand il y a des covariances)
phia	testInteractions	Permet de calculer tout type de contrastes, y compris les effets simples et les contrastes d’interaction
psych	describeBy	Obtenir les statistiques descriptives
psychometric	ICC1.lme	Permet de calculer l’ICC pour un modèle linéaire mixte
readxl	read_excel	Permet d’importer un fichier excel
sjstats	icc	Permet de calculer l’ICC pour un modèle linéaire mixte à partir d’un modèle ajusté avec le package lme4

|| || || ||

Solutions des exercices

Exercice 2

## Chargement des packages
library(nlme)
library(lme4)
library(lmerTest)
library(MuMIn)
library(HLMdiag)
library(phia)
library(psych)
library(reshape2) 
library(lattice)
library(HLMdiag)
library(nortest)
library(outliers)
library(psychometric)
library(sjstats)
library(emmeans)

Exercice 3

rat$Litter<-factor(rat$Litter) # le symbole dollar permet d'indiquer la variable "id" est dans le jeu de données "rat"
rat$Treatment<-factor(rat$Treatment)
rat$sex<-factor(rat$sex)

Exercice 4

rat2<-rat2[which(rat$Litter!=9),] 
nrow(rat2)

## [1] 305

Exercice 5

# distribution des résidus sur les données nettoyées avec lme
p1<-ggplot(rat.clean, aes(x=residus))+geom_histogram(aes(y=..density..))
p1<-p1+ stat_function(fun = dnorm, colour = "red",
                      args = list(mean = mean(rat.clean$residus, na.rm = TRUE),
                                  sd = sd(rat.clean$residus, na.rm = TRUE)))
p1<-p1+theme(plot.title = element_text(size = 12))+labs(x = "Distribution du résidu")
print(p1)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# distribution des résidus sur les données nettoyées avec lmer
p1<-ggplot(rat.clean2, aes(x=residus2))+geom_histogram(aes(y=..density..))
p1<-p1+ stat_function(fun = dnorm, colour = "red",
                      args = list(mean = mean(rat.clean2$residus, na.rm = TRUE),
                                  sd = sd(rat.clean2$residus, na.rm = TRUE)))
p1<-p1+theme(plot.title = element_text(size = 12))+labs(x = "Distribution du résidu")
print(p1)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Exercice 6

p2<-ggplot(rat.clean, aes(sample=residus))+stat_qq() 
p2<-p2+theme(plot.title = element_text(size = 12))+ggtitle("QQplot")
print(p2)

p2<-ggplot(rat.clean2, aes(sample=residus2))+stat_qq() 
p2<-p2+theme(plot.title = element_text(size = 12))+ggtitle("QQplot")
print(p2)

## Exercice 7

exo7<-lmer(weight~1+(1|Litter), data=rat.clean2)
pr01 <- profile(exo7)
xyplot(pr01, aspect = 1.3, layout=c(3,1)) # lignes verticales = IC à 0.5, 0.8, 0.9, 0.95, 0.99

xyplot(pr01, aspect = 1.3, layout=c(3,1), absVal=T)

Exercice 8

Modele.GLS0 <- gls(weight ~ Treatment + sex + Lsize +Treatment:sex, data = rat, method = "REML")
Modele.lme0<- lme(weight ~ Treatment + sex + Lsize +Treatment:sex,random=~1|Litter, data = rat, method = "REML")
anova(Modele.GLS0, Modele.lme0)

##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## Modele.GLS0     1  8 506.5099 536.5305 -245.2550                        
## Modele.lme0     2  9 419.1043 452.8775 -200.5522 1 vs 2 89.40562  <.0001

Exercice 9

Modele.lmer0<-lmer(weight~1+(1|Litter), data=rat)
vc<-VarCorr(Modele.lmer0, comp="Variance")
print(vc,comp=c("Variance","Std.Dev."))

##  Groups   Name        Variance Std.Dev.
##  Litter   (Intercept) 0.30037  0.54806 
##  Residual             0.19631  0.44307

Exercice 10

Modele.lme1<- lme(weight ~  1 , random=~1 | Litter,data = rat, method = "ML")
Modele.lme2<- lme(weight ~  sex + Lsize , random=~1 | Litter,data = rat, method = "ML")
Modele.lme3 <- lme(weight ~  sex + Lsize +Treatment , random=~1 | Litter,data = rat, method = "ML")
Modele.lme4 <- lme(weight ~  sex + Lsize +Treatment +Treatment:sex, random=~1 | Litter,data = rat, method = "ML")

anova(Modele.lme1,Modele.lme2,Modele.lme3,Modele.lme4)

##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## Modele.lme1     1  3 472.1246 483.4483 -233.0623                        
## Modele.lme2     2  5 407.4812 426.3540 -198.7406 1 vs 2 68.64344  <.0001
## Modele.lme3     3  7 392.7857 419.2076 -189.3929 2 vs 3 18.69549  0.0001
## Modele.lme4     4  9 395.8136 429.7845 -188.9068 3 vs 4  0.97212  0.6150

Exercice 11

r.squaredGLMM(Modele.lmer3)

##          R2m       R2c
## [1,] 0.43461 0.6238054

r.squaredLR(Modele.lmer3)

## [1] 0.5456628
## attr(,"adj.r.squared")
## [1] 0.6345789

Exercice 12

 contrasts(rat$Treatment)<-apply(contr.helmert(3), 2, rev)
 contrasts(rat$Treatment)

##         [,1] [,2]
## Control    0    2
## High       1   -1
## Low       -1   -1

 Modele.lmer.contraste <- lmer(weight ~  sex + Lsize +Treatment +Treatment:sex+ (1 | Litter),data = rat, REML = F)
 summary(Modele.lmer.contraste)

## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: weight ~ sex + Lsize + Treatment + Treatment:sex + (1 | Litter)
##    Data: rat
## 
##      AIC      BIC   logLik deviance df.resid 
##    395.8    429.8   -188.9    377.8      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.5188 -0.5045  0.0254  0.5880  3.0233 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Litter   (Intercept) 0.0807   0.2841  
##  Residual             0.1617   0.4021  
## Number of obs: 322, groups:  Litter, 27
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          7.51618    0.22025  41.63312  34.125  < 2e-16 ***
## sexMale              0.34584    0.05042 307.86545   6.859 3.81e-11 ***
## Lsize               -0.12821    0.01755  38.34969  -7.305 9.05e-09 ***
## Treatment1          -0.20814    0.08925  39.11922  -2.332 0.024935 *  
## Treatment2           0.19719    0.04682  38.13722   4.212 0.000149 ***
## sexMale:Treatment1  -0.01293    0.06623 310.89488  -0.195 0.845304    
## sexMale:Treatment2   0.03236    0.03281 304.43281   0.986 0.324756    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sexMal Lsize  Trtmn1 Trtmn2 sxM:T1
## sexMale     -0.071                                   
## Lsize       -0.954 -0.046                            
## Treatment1  -0.272 -0.133  0.342                     
## Treatment2   0.223  0.064 -0.265 -0.218              
## sxMl:Trtmn1 -0.044  0.304  0.008 -0.353  0.083       
## sxMl:Trtmn2 -0.016 -0.235  0.035  0.102 -0.386 -0.233

Exercice 13

 ez.reshape(data=ratbrain, 
            varying =list(c('Carb.BST','Carb.LS','Carb.VDB','Basal.BST','Basal.LS','Basal.VDB')), 
            v.names =c('Activite'),
            idvar =c('animal','IDeasy'),
            IV.names=list('Traitement','Region'), 
            IV.levels=list( c('carbachol','baseline'), c('BST','LS','VDB')))

## ez.reshape(data=ratbrain, varying =list(c('Carb.BST','Carb.LS','Carb.VDB','Basal.BST','Basal.LS','Basal.VDB')), v.names =c('Activite'),idvar =c('animal','IDeasy'),IV.names=list('Traitement','Region'), IV.levels=list( c('carbachol','baseline'), c('BST','LS','VDB')))

##     animal IDeasy      time Activite Traitement Region
## 1  R111097     p1  Carb.BST   371.71  carbachol    BST
## 2  R111397     p2  Carb.BST   492.58  carbachol    BST
## 3  R100797     p3  Carb.BST   664.72  carbachol    BST
## 4  R100997     p4  Carb.BST   515.29  carbachol    BST
## 5  R110597     p5  Carb.BST   589.25  carbachol    BST
## 6  R111097     p1   Carb.LS   302.02  carbachol     LS
## 7  R111397     p2   Carb.LS   355.74  carbachol     LS
## 8  R100797     p3   Carb.LS   587.10  carbachol     LS
## 9  R100997     p4   Carb.LS   437.56  carbachol     LS
## 10 R110597     p5   Carb.LS   493.93  carbachol     LS
## 11 R111097     p1  Carb.VDB   449.70  carbachol    VDB
## 12 R111397     p2  Carb.VDB   459.58  carbachol    VDB
## 13 R100797     p3  Carb.VDB   726.96  carbachol    VDB
## 14 R100997     p4  Carb.VDB   604.29  carbachol    VDB
## 15 R110597     p5  Carb.VDB   621.07  carbachol    VDB
## 16 R111097     p1 Basal.BST   366.19   baseline    BST
## 17 R111397     p2 Basal.BST   375.58   baseline    BST
## 18 R100797     p3 Basal.BST   458.16   baseline    BST
## 19 R100997     p4 Basal.BST   479.81   baseline    BST
## 20 R110597     p5 Basal.BST   462.79   baseline    BST
## 21 R111097     p1  Basal.LS   199.31   baseline     LS
## 22 R111397     p2  Basal.LS   204.85   baseline     LS
## 23 R100797     p3  Basal.LS   245.04   baseline     LS
## 24 R100997     p4  Basal.LS   261.19   baseline     LS
## 25 R110597     p5  Basal.LS   278.33   baseline     LS
## 26 R111097     p1 Basal.VDB   187.11   baseline    VDB
## 27 R111397     p2 Basal.VDB   179.38   baseline    VDB
## 28 R100797     p3 Basal.VDB   237.42   baseline    VDB
## 29 R100997     p4 Basal.VDB   195.51   baseline    VDB
## 30 R110597     p5 Basal.VDB   262.05   baseline    VDB

Exercice 14

 modele0<-lme(Activite~Traitement*Region, random=~1|animal, data= ratbrain.long, method="ML")

Exercice 15

modele4a<-lme(Activite~Traitement*Region, 
        random=~1|animal/Traitement, data= ratbrain.long, method="REML")
modele4b<-lme(Activite~Traitement*Region, 
        random=~1|animal/Traitement, data= ratbrain.long, method="REML",
        weights = varIdent(form = ~ 1 | Traitement))

 anova(modele4a, modele4b)

##          Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## modele4a     1  9 274.7148 285.3173 -128.3574                         
## modele4b     2 10 276.5303 288.3109 -128.2652 1 vs 2 0.1844388  0.6676

Exercice 16

Il existe plusieurs manières d’obtenir les contrastes. La manière la plus simple en l’occurrence est d’utiliser la fonction summary.

summary(modele1b)

## Linear mixed-effects model fit by maximum likelihood
##  Data: ratbrain.long 
##       AIC      BIC    logLik
##   312.723 326.7349 -146.3615
## 
## Random effects:
##  Formula: ~Traitement | animal
##  Structure: General positive-definite, Log-Cholesky parametrization
##                    StdDev   Corr  
## (Intercept)        98.94686 (Intr)
## Traitementbaseline 71.39379 -0.981
## Residual           20.76334       
## 
## Fixed effects: Activite ~ Traitement * Region 
##                                 Value Std.Error DF    t-value p-value
## (Intercept)                   526.710  50.55096 20  10.419387  0.0000
## Traitementbaseline            -98.204  38.59827 20  -2.544259  0.0193
## RegionLS                      -91.440  14.68190 20  -6.228076  0.0000
## RegionVDB                      45.610  14.68190 20   3.106546  0.0056
## Traitementbaseline:RegionLS   -99.322  20.76334 20  -4.783526  0.0001
## Traitementbaseline:RegionVDB -261.822  20.76334 20 -12.609819  0.0000
##  Correlation: 
##                              (Intr) Trtmnt RegnLS RgnVDB Tr:RLS
## Traitementbaseline           -0.943                            
## RegionLS                     -0.145  0.190                     
## RegionVDB                    -0.145  0.190  0.500              
## Traitementbaseline:RegionLS   0.103 -0.269 -0.707 -0.354       
## Traitementbaseline:RegionVDB  0.103 -0.269 -0.354 -0.707  0.500
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -1.745380834 -0.451592586 -0.007145421  0.518925580  2.040618276 
## 
## Number of Observations: 30
## Number of Groups: 5

emmeans.out<- emmeans(modele1b, ~ Traitement*Region, weights="show.levels")
plot(emmeans.out)

emmeans.out<-summary(emmeans.out)
interaction.plot(ratbrain.long$Region, ratbrain.long$Traitement, ratbrain.long$Activite, las=1)

pd <- position_dodge(0.1) 
ggplot(emmeans.out, aes(x=Region, y=emmean, colour=Traitement, group=Traitement)) + 
    geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, size=3)

Exercice 17

autisme$lang<-ifelse(autisme$sicdegp==1, "faible", ifelse(autisme$sicdegp==2, "moyen", "eleve"))
autisme$lang<-factor(autisme$lang)
autisme$age.f<-factor(autisme$age)

Exercice 18

interaction.plot(autisme$age, autisme$childid, autisme$vsae, las=1)

Exercice 19 :

autisme2a<-update(autisme2, random=~age.2-1)
anova(autisme2, autisme2a)

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme2      1 13 4641.276 4698.457 -2307.638                        
## autisme2a     2 11 4721.203 4769.587 -2349.601 1 vs 2 83.92683  <.0001

Exercice 20

## vsae ~ age.2 + I(age.2^2) + lang + age.2:lang + I(age.2^2):lang
## <environment: 0x000000003b6bbc20>

## vsae ~ age.2 + I(age.2^2) + lang + age.2:lang
## <environment: 0x000000003e1b8f78>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3      1 13 4636.444 4693.819 -2305.222                        
## autisme3a     2 11 4634.314 4682.862 -2306.157 1 vs 2 1.869704  0.3926

## vsae ~ age.2 + I(age.2^2) + lang
## <environment: 0x000000003e1b8f78>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3b     2  9 4653.696 4693.417 -2317.848 1 vs 2 23.38232  <.0001

## vsae ~ age.2 + I(age.2^2) + age.2:lang
## <environment: 0x000000003e1b8f78>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3c     2  9 4653.929 4693.650 -2317.964 1 vs 2 23.61505  <.0001

## vsae ~ age.2 + lang + age.2:lang
## <environment: 0x000000003e1b8f78>

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## autisme3a     1 11 4634.314 4682.862 -2306.157                        
## autisme3d     2 10 4638.766 4682.900 -2309.383 1 vs 2 6.451775  0.0111

## vsae ~ I(age.2^2) + lang + lang:age.2
## <environment: 0x000000003e1b8f78>

##           Model df      AIC      BIC    logLik
## autisme3a     1 11 4634.314 4682.862 -2306.157
## autisme3e     2 11 4634.314 4682.862 -2306.157

##                 numDF denDF  F-value p-value
## (Intercept)         1   446 949.1017  <.0001
## age.2               1   446 145.9686  <.0001
## I(age.2^2)          1   446   6.1842  0.0133
## lang                2   155  23.0890  <.0001
## age.2:lang          2   446  12.6204  <.0001
## I(age.2^2):lang     2   446   0.9265  0.3967

Le modèle 3a est le modèle qu’on garde. Il faut constater que le fait d’avoir ou on l’âge ne change pas la vraisemblance ici (sans que l’explication ne soit claire à mes yeux), mais que l’anova permet de voir que l’âge a un effet signficatif.

Supplément à la formation

Modèles linéaires mixtes et ressources de calcul

Les effets croisés requièrent plus de ressources de calcul que les effets emboîtés. Ceci s’explique par le fait que les matrices du design ne sont plus diagonales. Il est préférable dans ce cas d’utiliser la fonction lmer du package lme4 qui est optimsée pour ce type de design.

str(Dyestuff)

## 'data.frame':    30 obs. of  2 variables:
##  $ Batch: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ Yield: num  1545 1440 1440 1520 1580 ...

fm01 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff)
fm01@pp$LamtUt

## 6 x 30 sparse Matrix of class "dgCMatrix"

##    [[ suppressing 30 column names '1', '2', '3' ... ]]

##                                                                           
## [1,] 0.8483238 0.8483238 0.8483238 0.8483238 0.8483238 .         .        
## [2,] .         .         .         .         .         0.8483238 0.8483238
## [3,] .         .         .         .         .         .         .        
## [4,] .         .         .         .         .         .         .        
## [5,] .         .         .         .         .         .         .        
## [6,] .         .         .         .         .         .         .        
##                                                                           
## [1,] .         .         .         .         .         .         .        
## [2,] 0.8483238 0.8483238 0.8483238 .         .         .         .        
## [3,] .         .         .         0.8483238 0.8483238 0.8483238 0.8483238
## [4,] .         .         .         .         .         .         .        
## [5,] .         .         .         .         .         .         .        
## [6,] .         .         .         .         .         .         .        
##                                                                           
## [1,] .         .         .         .         .         .         .        
## [2,] .         .         .         .         .         .         .        
## [3,] 0.8483238 .         .         .         .         .         .        
## [4,] .         0.8483238 0.8483238 0.8483238 0.8483238 0.8483238 .        
## [5,] .         .         .         .         .         .         0.8483238
## [6,] .         .         .         .         .         .         .        
##                                                                           
## [1,] .         .         .         .         .         .         .        
## [2,] .         .         .         .         .         .         .        
## [3,] .         .         .         .         .         .         .        
## [4,] .         .         .         .         .         .         .        
## [5,] 0.8483238 0.8483238 0.8483238 0.8483238 .         .         .        
## [6,] .         .         .         .         0.8483238 0.8483238 0.8483238
##                         
## [1,] .         .        
## [2,] .         .        
## [3,] .         .        
## [4,] .         .        
## [5,] .         .        
## [6,] 0.8483238 0.8483238

Dans l’exemple ci-dessus, on a plusieurs mesures pour le même individu, mais aucun facteur en mesure répétée n’a été introduit. On travaille donc sur la diagonale uniquement (la valeur est la même pour toutes les mesures).

Ce qu’il faut constater ici est que seule la variance est estimée ce qui permet au calcul d’être relativement rapide. En revanche, quand les éléments en dehors de la diagonale ne valent pas 0, le temps de calcul va être considérablement augmenter.

Effets aléatoires et coefficients de régression.

Bien que, pour les effets aléatoires, les coefficients de régressions ne soient pas explicitement estimés, il est possible d’en calculer un. C’est ce qu’on appelle le BLUP (best linear unbiased predictors). L’intérêt est de pouvoir faire des inférences sur la variabilité des effets aléatoires.

Bates se montre néanmoins critique sur l’intérêt de cet indicateur.

Les valeurs prédites : valeurs prédites conditonnelles et inconditonnelles.

La prédiction est une fonction essentielle des statistiques, en particulier des modèles linéaires. La prédiction consiste à déterminer la valeur qu’aurait dû obtenir une observation particulière. Ces valeurs prédites sont obtenues par la droite de régression.

Dans le cas du modèle linéaire mixte, il est possible de prédire la valeur que devrait avoir chaque observation de deux manières: soit on l’estime à l’intérieur de chacune des modalités de la variable aléatoire, sont on l’estime pour tous les groupes confondus. Dans le premier cas, on parle de valeurs prédites conditionnelles et, dans le second de valeurs prédites inconditionnelles.

Dans le premier cas, la formule prend la forme suivante :

\[ Poids = 7.723941 + 0.17 \times SEXE + (-0.13) \times Taille.Portée + 0.43 \times Traitement_{haut} + (-0.001) \times Traitement_{bas} + û_j \]

Dans le second cas, la formule prend la forme suivante :

\[ Poids = 7.723941 + 0.17 \times SEXE + (-0.13) \times Taille.Portée + 0.43 \times Traitement_{haut} + (-0.001) \times Traitement_{bas} \]

Ainsi, la différence entre les deux est la prise en compte ou non de l’effet aléatoire.

Pour le formuler autrement, il va y avoir un coefficient par portée. Comme nous avons 27 portées, nous aurons 27 coefficients. Ces coefficients s’appellent les meilleurs prédicteurs linéaires non biaisés (BLUP).

On les obtient avec le packages nlme par la fonction “random.effects”

 random.effects(Modele.lme)

##    (Intercept)
## 1   0.17805612
## 2  -0.07071525
## 3  -0.18289450
## 4  -0.05418319
## 5   0.34615302
## 6  -0.05363520
## 7   0.39046482
## 8  -0.02968245
## 9  -0.60785892
## 10  0.08429556
## 11  0.03832949
## 12  0.02490328
## 13 -0.36672606
## 14  0.01557233
## 15  0.03327425
## 16  0.07051165
## 17 -0.39760275
## 18  0.43692016
## 19 -0.18154659
## 20  0.32636423
## 21 -0.29062992
## 22 -0.49231419
## 23  0.26451251
## 24  0.18828098
## 25  0.23099259
## 26  0.24197633
## 27 -0.14281830

et par la fonction “ranef” dans le package lme4.

 ranef(Modele.lmer)

## $Litter
##    (Intercept)
## 1   0.17480024
## 2  -0.07362296
## 3  -0.17490203
## 4  -0.05376249
## 5   0.34446954
## 6  -0.05480208
## 7   0.39153639
## 8  -0.02616704
## 9  -0.61772107
## 10  0.09017150
## 11  0.04136696
## 12  0.02072931
## 13 -0.35981737
## 14  0.01847368
## 15  0.03549783
## 16  0.06416884
## 17 -0.39636862
## 18  0.42802096
## 19 -0.18110866
## 20  0.32903708
## 21 -0.27813902
## 22 -0.49096621
## 23  0.26053477
## 24  0.17537803
## 25  0.23748828
## 26  0.22966912
## 27 -0.13396497
## 
## with conditional variances for "Litter"

Ainsi, une différence majeure entre le MLM et le modèle linéaire est que, contrairement au modèle linéaire, on peut prédire les valeurs de l’effet aléatoire. Il s’agit non pas de valeurs fixes comme les bêtas mais de valeurs aléatoires qui suivent une distribution multinormale. L’objectif est donc de prédire les valeurs des effets aléatoires plutôt que de les estimer.

Contrairement aux effets fixes, nous ne voulons pas prédire la moyenne (la valeur attendue) mais on peut utiliser l’ensemble des informations pour prédire les valeurs attendues des effets aléatoires.

Evidemment, il n’est pas nécessaire de calculer ces valeurs à la main. On peut les obtenir en allant chercher dans le modèle les valeurs ajustées.

 Modele.lme$fitted[1:10,]

##       fixed   Litter
## 1  6.761836 6.939893
## 2  6.761836 6.939893
## 3  6.761836 6.939893
## 4  6.761836 6.939893
## 5  6.761836 6.939893
## 6  6.761836 6.939893
## 7  6.761836 6.939893
## 8  6.761836 6.939893
## 9  6.402755 6.580811
## 10 6.402755 6.580811

Dans notre exemple, voici les valeurs pour le rat numéro 1.

 rat[1,]

## # A tibble: 1 x 8
##   id    weight sex   Litter Lsize Treatment residus residus2
##   <fct>  <dbl> <fct> <fct>  <dbl> <fct>       <dbl>    <dbl>
## 1 1        6.6 Male  1         12 Control    -0.670   -0.884

Les valeurs prédites pour ce rat sont donc :

SEXE1<-2
Taille<-12
Haut<-0
Bas<-0
Liter<-1
7.723941 + 0.1717157*SEXE1 + (-0.1306805*(Taille-2)) +0.43*Haut+(-0.001)*Bas # valeur prédite inconditionnelle

## [1] 6.760567

7.887512 + 0.1717157*SEXE1 + (-0.1306805*(Taille-2)) +0.43*Haut+(-0.001)*Bas # valeur prédite conditionnelle

## [1] 6.924138

Remarquez que l’intercept a changé dans le modèle conditionnel. Il s’agit de l’intercept spécifique à la portée. Ce coefficient peut être obtenu de la manière suivante :

coef(Modele.lme)

##    (Intercept) Treatmentcont1 Treatmentcont2   sexMale      Lsize
## 1     7.699782      0.2148158      0.2142509 0.3590819 -0.1290031
## 2     7.451010      0.2148158      0.2142509 0.3590819 -0.1290031
## 3     7.338831      0.2148158      0.2142509 0.3590819 -0.1290031
## 4     7.467542      0.2148158      0.2142509 0.3590819 -0.1290031
## 5     7.867878      0.2148158      0.2142509 0.3590819 -0.1290031
## 6     7.468090      0.2148158      0.2142509 0.3590819 -0.1290031
## 7     7.912190      0.2148158      0.2142509 0.3590819 -0.1290031
## 8     7.492043      0.2148158      0.2142509 0.3590819 -0.1290031
## 9     6.913867      0.2148158      0.2142509 0.3590819 -0.1290031
## 10    7.606021      0.2148158      0.2142509 0.3590819 -0.1290031
## 11    7.560055      0.2148158      0.2142509 0.3590819 -0.1290031
## 12    7.546629      0.2148158      0.2142509 0.3590819 -0.1290031
## 13    7.154999      0.2148158      0.2142509 0.3590819 -0.1290031
## 14    7.537298      0.2148158      0.2142509 0.3590819 -0.1290031
## 15    7.555000      0.2148158      0.2142509 0.3590819 -0.1290031
## 16    7.592237      0.2148158      0.2142509 0.3590819 -0.1290031
## 17    7.124123      0.2148158      0.2142509 0.3590819 -0.1290031
## 18    7.958646      0.2148158      0.2142509 0.3590819 -0.1290031
## 19    7.340179      0.2148158      0.2142509 0.3590819 -0.1290031
## 20    7.848090      0.2148158      0.2142509 0.3590819 -0.1290031
## 21    7.231096      0.2148158      0.2142509 0.3590819 -0.1290031
## 22    7.029411      0.2148158      0.2142509 0.3590819 -0.1290031
## 23    7.786238      0.2148158      0.2142509 0.3590819 -0.1290031
## 24    7.710006      0.2148158      0.2142509 0.3590819 -0.1290031
## 25    7.752718      0.2148158      0.2142509 0.3590819 -0.1290031
## 26    7.763702      0.2148158      0.2142509 0.3590819 -0.1290031
## 27    7.378907      0.2148158      0.2142509 0.3590819 -0.1290031

Le maximum de vraisemblance.

En général, le maximum de vraisemblance est un méthode pour obtenir les estimations de paramètres inconnus en optimisant la fonction de vraisemblance.

Le principe la vraisemblance est d’estimer la probabilité d’observer les données de manière itérative. Ainsi, on part d’une valeur “aléatoire” et on va déplacer cette valeur aléatoire de sorte à augmenter la probabilité d’observer cette distribution de données par simulation. Le maximum de vraisemblance est la valeur de la probabilité pour laquelle les paramètre seront optimisés pour avoir la probabilité maximale. Quand cette probabilité a atteint un maximum (moyennant une certaine tolérance), on dit que le modèle a convergé.

Pour une variable numérique, la vraisemblance est maximale lorsque la moyenne et l’écart-type sont identiques à la moyenne et l’écart-type calculés.

mean(rat$weight)

## [1] 6.080963

sd(rat$weight)

## [1] 0.6474272

Ainsi, pour chaque observation d’une variable numérique, on peut estimer la densité de probabilité de cette observation sur la base d’une distribution normale (pour simplifier un peu, la hauteur de la courbe normale).

rat$weight[1:10]

##  [1] 6.60 7.40 7.15 7.24 7.10 6.04 6.98 7.05 6.95 6.29

dnorm(rat$weight, mean= 6.10 , sd= 0.64, log=F)[1:10]

##  [1] 0.45940311 0.07921103 0.16227276 0.12757164 0.18390158 0.62061400
##  [7] 0.24220666 0.20714505 0.25804806 0.59647448

La fonction de vraisemblance est obtenue en appliquant le logarithme à chacune des observations.

dnorm(rat$weight, mean= 6.10 , sd= 0.64, log=T)[1:10]

##  [1] -0.7778272 -2.5356397 -1.8184766 -2.0590772 -1.6933546 -0.4770460
##  [7] -1.4179639 -1.5743360 -1.3546094 -0.5167188

On obtient la vraisemblance en additionnant le logarithme des probabilités d’observer cette distribution.

sum(dnorm(rat$weight, mean= 6.10 , sd= 0.64, log=T))

## [1] -316.583

Cette valeur se distribue en suivant une distribution de \(\chi^2\).

Dans le cas du modèle linéaire mixte, la fonction de la densité de probabilité pour une distribution normale multivariée est obtenue par la formule suivante :

\[ f(Y_i|\beta, \theta) = (2\pi)^{-n_i/2} det(V_i)^{-1/2} exp(-0.5 \times (Y_i-X_i \beta)'V_i^{-1}(Y_i-X_i \beta)) \] Sur la base de cette fonction de densité de probabilité, on peut déterminer la contribution de la fonction de vraisemblance pour un participant i donné.

\[ L_i(Y_i|\beta, \theta;y_i) = (2\pi)^{-n_i/2} det(V_i)^{-1/2} exp(-0.5 \times (Y_i-X_i \beta)'V_i^{-1}(Y_i-X_i \beta)) \]

Et la fonction de vraisemblance est le produit des m contributions indépendantes des individus (i = 1…m).

\[ L_i(Y_i|\beta, \theta) =\prod (2\pi)^{-n_i/2} det(V_i)^{-1/2} exp(-0.5 \times (Y_i-X_i \beta)'V_i^{-1}(Y_i-X_i \beta)) \]

Ce qu’il faut identifier dans la formule, c’est que cette estimation dépend de la matrice de Variance/covariance et de l’optimisation des moindres carrés.

modelML <- lme(weight ~ Treatment + sex + Lsize +Treatment:sex, random = ~ 1 | Litter,data = rat, method = "ML")
 summary(modelML)$logLik

## [1] -188.9068

Dans ce modèle, la vraisemblance est égale à -188.9068.

Pour obtenir la déviance, on multiplie cette vraisemblance par -2.

-2*-188.9068

## [1] 377.8136

Un des biais du maximum de vraisemblance est que les estimations de certains paramètres pour les effets aléatoires (theta) ne prennent pas en compte la perte de degrés de liberté qui résulte de l’estimation des paramètres des effets fixes. Un moyen d’éliminer ce biais est d’utiliser le maximum de vraisemblance restreint (REML)

Cette méthode a été introduite pour estimer des composantes de variance/covariance dans le contexte de plan qui ne sont pas équilibrés, ou incomplets. On préfère généralement le REML parce qu’il produit des estimations non biaisées des paramètres de covariances en prenant en compte la perte de ddl qui résulte de l’estimation des effets fixes.

Pour ces raisons, cet estimateur est souvent préféré au maximum de vraisemblance. Cependant, Bates souligne que la distinction entre les deux est probablement moins importante que ce que la plupart des personnes l’imagine (p. 110).

modelREML <- lme(weight ~ Treatment + sex + Lsize +Treatment:sex, random = ~ 1 | Litter,data = rat, method = "REML")
 summary(modelREML)$logLik

## [1] -204.1357

Dans ce modèle, le log de vraisemblance restreint est égal à -204.8288

Le ratio de vraisemblance (LRT)

Le LRT nécessite que le modèle emboîté et le modèle de référence soient ajustés sur le même jeu de données. La différence de vraisemblance entre deux modèles se distribue approximative comme un \(\chi^2\).

Par exemple, nous pouvons savoir si notre modèle dans son ensemble est plus vraisemblable que la modélisation de la constante uniquement.

Nous devons donc modéliser la constante uniquement. Pour cela, l’équation prend la forme :

\[ VD \sim 1 \]

Modele.LRT0 <- lme(weight ~ 1, random = ~ 1 | Litter,data = rat, method = "ML")

Nous devons avoir le modèle que nous voulons tester.

Modele.LRT.complet <- lme(weight ~ Treatment + sex + Lsize +
                      Treatment:sex, random = ~ 1 | Litter,data = rat, method = "ML")

Nous pouvons à présent comparer notre modèle par rapport à un modèle où seul la constante est modélisée en utilisant la fonction anova.

results<-anova(Modele.LRT0, Modele.LRT.complet)
results

##                    Model df      AIC      BIC    logLik   Test  L.Ratio
## Modele.LRT0            1  3 472.1246 483.4483 -233.0623                
## Modele.LRT.complet     2  9 395.8136 429.7845 -188.9068 1 vs 2 88.31105
##                    p-value
## Modele.LRT0               
## Modele.LRT.complet  <.0001

Dans les résultats, on retrouve le log de vraisemblance du modèle complet vaut \(-188.9067913\).De la même manière la vraisemblance du modèle où seul l’intercept est modélisé vaut \(-233.0623144\).

La significativité est obtenue en calculant la déviance, c’est-à-dire en multipliant la différence entre ces deux logarithmes de vraisemblance par -2. La déviance suit une distribution \(\chi^2\) ayant comme degrés de liberté le nombre de paramètres estimés.

-2 * (-233.0623-(-188.9068))

## [1] 88.311

Ainsi, un \(\chi^2\) de 88.311 ayant 6 (9-3) degrés de liberté est significatif à un seuil inférieur à 0.0001.

Dans cet exemple, les 6 paramètres estimés sont : - Le bêta pour l’effet du sexe. - Le bêta pour l’effet de la taille de la portée. - Deux bêtas pour l’effet du Traitement. - Deux bêtas pour l’effet d’interaction entre le sexe et le traitement.

Selon West, Welch et Gatecki (2015), il est préférable d’utiliser le maximum de vraisemblance pour les effets fixes et le maximum de vraisemblance restreint pour les effets aléatoires, en particulier si la taille de l’échantillon est faible car le REML réduit les biais inhérent à l’estimation du ML pour les paramètres de covariances.

Ainsi, on peut utiliser par exemple ce ratio de vraisemblance sur les effets aléaoires pour :

• tester sur les variances résiduelles sont hétérogènes vs homogènes.

• tester si les covariances entre deux effets aléatoires sont nulles.

• savoir si cela a du sens de garder tous les effets aléatoires. En fait, on ne va pas tester directement les effets aléatoire, mais on teste si les covariance des effets aléatoires sont égales à 0. S’il n’y a qu’un seul effet aléatoire, on peut se questionner sur la possibilité d’omettre cet effet aléatoire.

On obtient les mêmes résultats avec le package lme4.

 model3.1.lmer0 <- lmer(weight ~ 1 + (1 | Litter),data = rat, REML=F)
 model3.1.lmer <- lmer(weight ~ Treatment + sex + Lsize +Treatment:sex +( 1 | Litter),data = rat, REML=F)
 anova( model3.1.lmer0,  model3.1.lmer)

## Data: rat
## Models:
## model3.1.lmer0: weight ~ 1 + (1 | Litter)
## model3.1.lmer: weight ~ Treatment + sex + Lsize + Treatment:sex + (1 | Litter)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model3.1.lmer0  3 472.12 483.45 -233.06   466.12                         
## model3.1.lmer   9 395.81 429.78 -188.91   377.81 88.311      6  < 2.2e-16
##                   
## model3.1.lmer0    
## model3.1.lmer  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Les critères d’information.

Les critères d’information fournissent une manière d’évaluer l’ajustement d’un modèle sur la base la valeur optimale du log de vraisemblance tout en pénalisant les modèles qui ne sont pas parcimonieux.

L’intérêt des critères d’information est qu’ils permettent de comparer n’importe quels modèles ajustés sur le même jeu de données. Il n’est donc pas nécessaire que les modèles soient emboîtés.

Il existe deux critères d’information qui sont particulièrement utilisés : le AIC (Critère d’Informations d’Akaike) et le BIC (Critère d’Information de Bayes).

Ils sont fournis quand on demande le résumé d’un modèle.

summary(model3.1.lmer)$AIC

##       AIC       BIC    logLik  deviance  df.resid 
##  395.8136  429.7845 -188.9068  377.8136  313.0000

Dans les deux cas, la logique est la même : on va s’appuyer sur le calcul de la vraisemblance, mais en pénalisant le modèle à chaque paramètre estimé.

Pour le AIC, la déviance va être pénalisée en additionnant 2 fois le nombre de prédicteur.

\[ AIC=-2×l(β ̂,θ ̂ )+2p \]

Dans notre modèle, le log de vraisemblance vaut -188.9068, la déviance vaut donc -2*-188.9068, ce qui donne 377.8136. Les degrés de liberté résiduels sont de 313 alors que le nombre d’observations est de 322. La différence entre 322 et 313 est de 9. Il y a donc 9 paramètres qui ont été estimés. Le AIC vaut donc :

-2*summary(model3.1.lmer)$logLik+2*9

## 'log Lik.' 395.8136 (df=9)

Pour le critère d’information de Bayes, le principe est assez similaire, la différence porte sur la manière de calculer la pénalité pour le manque de parcimonie. Ainsi, le nombre de paramètres estimés sera multiplié par le logarithme de la taille de l’échantillon.

\[BIC=-2×l(β ̂,θ ̂ )+p×lna(n) \]

-2*summary(model3.1.lmer)$logLik+9*log(322, base = exp(1))

## 'log Lik.' 429.7845 (df=9)

Puisque ces deux critères d’informations s’appuient sur la déviance pour le calcul, à laquelle on va ajouter une pénalité aux modèles pour lesquels de nombreux paramètres sont estimés, cela signifie que, comme pour la déviance où plus la valeur est faible, plus les données sont ajustées, plus la valeur pour ces deux critères d’informations est petite, meilleur est l’ajustement. Dans notre exemple, tant le AIC que le BIC sont meilleurs pour le modèle ajusté par rapport au modèle où seul l’intercept est modélisé.

 anova( model3.1.lmer0,  model3.1.lmer)

## Data: rat
## Models:
## model3.1.lmer0: weight ~ 1 + (1 | Litter)
## model3.1.lmer: weight ~ Treatment + sex + Lsize + Treatment:sex + (1 | Litter)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## model3.1.lmer0  3 472.12 483.45 -233.06   466.12                         
## model3.1.lmer   9 395.81 429.78 -188.91   377.81 88.311      6  < 2.2e-16
##                   
## model3.1.lmer0    
## model3.1.lmer  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pour terminer avec ces critères d’informations, il est important de comprendre que, de manière théorique, le AIC et le BIC peuvent être calculés sur la base du ML ou du REML d’un modèle ajusté. Cependant, dans le package ‘lme4’, ils ne sont pas calculés avec le REML étant donné qu’il ne s’agit pas d’un vrai maximum de vraisemblance.

Les fonctions d’optimisations.

Comme le maximum de vraisemblance est un processus itératif, son estimation est intensive en termes de ressources de calcul. Si on prenait des valeurs au hasard, l’estimation du maximum de vraisemblance ne serait pas optimale. Il est donc nécessaire d’optimiser le calcul de ces fonctions.

Il existe 3 grands algorithmes d’optimisation : la maximisation attendue (expectation-maximization - EM), l’algorithme de Newton-Raphson (N-R) et les scores de fisher (Fisher scoring).

• L’algorithme EM : idéal pour trouver de bonnes valeurs de départ et pour maximiser des fonctions de vraisemblance qui sont compliquées. Pour converger, EM a besoin d’un jeu de données équilibré (i.e., avoir le même nombre d’observations dans chaque condition). Si ce n’est pas le cas, il va augmenter le jeu de données de manière artificielle en utilisant la somme des sommes des carrés. L’a raison est’idée sous-jacente de cet algorithme est qu’il est plus simple d’optimiser la vraisemblance sur la base de données complètes que sur les données observées. Le problème est que cet algorithme est lent à converger et est généralement trop optimiste puisqu’on utilise les données complètes plutôt qu’observées. On utilise rarement cet algorithme pour l’estimation de la vraisembance, il est en revanche utile pour optimiser les valeurs de départ.

• L’algorithme de N-R est l’algorithme le plus souvent utilisé. Cet algorithme minimise \(-2 * log (vraisemblance)\). Si cette méthode consomme plus de temps à chaque itération, on atteint la convergence habituellement avec moins d’itérations.

• La méthode des scores de Fisher se distingue de la méthode précédente en utilisant la matrice attendue plutôt qu’observée. Cette méthode est généralement plus stable et va probablement converger mais elle n’est pas recommandée pour obtenir les estimateurs. En effet, la difficulté est de déterminer les valeurs attendues. Donc, pour éviter ces problèmes, il vaut mieux utiliser la méthode N-R

R va utiliser l’algorithme EM pour fixer les valeurs de départ et l’algorithme NR pour la convergence.

La raison pour laquelle c’est important de comprendre à quoi cela correspond est qu’il arrive parfois que le modèle ne converge pas et il est possible de paramétrer ces fonctions d’optimisation pour aider le modèle à converger.

 lmeControl() # fonction permettant de paramétrer les fonctions d'optimisation. 

## $maxIter
## [1] 50
## 
## $msMaxIter
## [1] 50
## 
## $tolerance
## [1] 1e-06
## 
## $niterEM
## [1] 25
## 
## $msMaxEval
## [1] 200
## 
## $msTol
## [1] 1e-07
## 
## $msVerbose
## [1] FALSE
## 
## $returnObject
## [1] FALSE
## 
## $gradHess
## [1] TRUE
## 
## $apVar
## [1] TRUE
## 
## $.relStep
## [1] 6.055454e-06
## 
## $opt
## [1] "nlminb"
## 
## $optimMethod
## [1] "BFGS"
## 
## $minAbsParApVar
## [1] 0.05
## 
## $natural
## [1] TRUE
## 
## $sigma
## [1] 0
## 
## $allow.n.lt.q
## [1] FALSE

Par défaut, on identifie ici que le nombre maximum d’itération du modèle est de 50 et que l’estimations des valeurs de départ doivent s’obtenir avec un maximum de 25 itérations. Enfin, on tolère que la vraisemblance a atteint son maximum si entre deux itérations, l’amélioration de la vraisemblance est inférieure à \(10^{-6}\).

Il est possible de modifier ces paramètres de la manière suivante.

 lmeControl(maxIter=100, tolerance=10^-5) 

## $maxIter
## [1] 100
## 
## $msMaxIter
## [1] 50
## 
## $tolerance
## [1] 1e-05
## 
## $niterEM
## [1] 25
## 
## $msMaxEval
## [1] 200
## 
## $msTol
## [1] 1e-07
## 
## $msVerbose
## [1] FALSE
## 
## $returnObject
## [1] FALSE
## 
## $gradHess
## [1] TRUE
## 
## $apVar
## [1] TRUE
## 
## $.relStep
## [1] 6.055454e-06
## 
## $opt
## [1] "nlminb"
## 
## $optimMethod
## [1] "BFGS"
## 
## $minAbsParApVar
## [1] 0.05
## 
## $natural
## [1] TRUE
## 
## $sigma
## [1] 0
## 
## $allow.n.lt.q
## [1] FALSE

Si devoir utiliser plus d’itérations pour converger est tout à fait acceptable puisque le seul impact est simplement une augmentation du temps de calcul accordé à l’ordinateur pour trouver la fonction d’optimisation, diminuer la tolérance est plus délicat puisque cela diminue la précision. Cependant, à ma connaissance, il n’existe pas de règle empirique fixant les limites.

Problèmes computationnels dans l’estimation des paramètres de covariances.

Il arrive que, lorsque les corrélations entre les observations sont fortes, ou que le modèle ait été mal spécifié, l’estimation pour les paramètres theta des facteurs aléatoires peut se trouver en dehors des limites de l’espace du paramètre, ce qui a pour conséquence que la matrice n’est pas définie positive.

Lorsqu’on réalise un modèle linéaire mixte, il faut faire attention aux messages d’avertissement/erreur indiquant que la matrice n’est pas définie positive. Cela signifie que les valeurs calculées sont très fortement biaisées et ne devraient pas être interprétées.

Pour contourner ce problème :

• on peut choisir des valeurs de départ alternatives pour l’estimation des paramètres de covariance (le paramètre “start” de la fonction lmer)

• Changer l’échelle des données, surtout si les données ne sont pas distribuées en suivant la multinormalité. Cela survient aussi quand on a des variances énormes. Par exemple, on mesure un temps un minutes alors que la durée mesurée est de plusieurs jours. Le fait de passer en heures peut lever la difficulté de convergence.

• Simplifier le modèle en supprimant les effets aléatoires qui ne sont nécessaires. On peut aussi simplifier le modèle sur les effets fixes. Généralement, il est recommandé de supprimer les termes de niveaux supérieurs, comme les interactions de niveaux supérieurs. Dans tous les cas, il faut se questionner sur les effets qui doivent être maintenus et ceux qui ne sont pas nécessaires, car il faut pouvoir le justifier.

Il existe d’autres méthodes, mais peu importe la méthode, aucune ne permet de garantir que la modèle va converger.

Stratégie de construction de modèle

Dans un modèle linéaire mixte, la difficulté réside dans la capacité à choisir le modèle le plus simple qui s’ajuste le mieux aux données.

Le processus de construction de modèle est un processus itératif qui nécessite une série d’étapes dans l’ajustement du modèle. Peu importe la stratégie, il n’y en a pas qui puisse s’appliquer à toutes les sitations.

La stratégie Top-Down

1. On débute avec des moyennes bien définies pour le modèle (ie, un maximum d’effets fixes). C’est ce qu’on appelle un modèle avec un structure saturée.

2. On sélectionne une structure pour les effets aléatoires dans le modèle. Pour cela, il faut choisir les variables qui doivent être introduites comme variables aléatoires dans le modèle. Pour tester cela, on peut utiliser le ratio de vraisemblance avec comme alogirthme le maximum de vraisemblance restreint (REML).

3. Choisir la structure de covariance pour les résidus. A cette étape, il s’agit de déterminer la structure de la variance résiduelle (utile uniquement pour des mesures répétées/longitudinales).

4. Réduire le modèle. Il s’agit de vérifier en utilisant le test approprié si tous les paramètres des effets fixes sont nécessaires dans le modèle.

On répète les opérations 2 à 4 autant de fois que nécessaire pour obtenir le modèle final.

La stratégie pas-à-pas ascendant.

Les étapes à suivre dans la stratégie pas-à-pas ascendant sont les suivantes :

1. On commence avec un modèle inconditionnel. Ainsi, on ne modélise que l’intercept comme effet fixe et les effets aléatoires.

2. Construire le modèle en ajoutant la covariable de niveau 1 au premier niveau du modèle. Au niveau 2 du modèle, il faut se questionner sur l’intérêt d’ajouter les effets aléatoires pour les coefficients de niveau 1.

3. On construit ensuite le modèle en ajoutant les covariables de niveau 2. Pour un modèle à 3 niveaux, il faut se questionner sur l’intérêt d’ajouter les effets aléatoires pour les coefficients de niveau 2. A cette étape, Les covariables de niveau 2 et les effets fixes qui y sont associés peuvent être ajouter au niveau 2 du modèle. Ces covariables de niveau 2 pourraient expliquer certains effets aléatoires du niveau 1 qui seraient ainsi pris en compte au niveau des effets aléatoires de niveau 2.

4. l’opération est répétée au niveau 3 pour un modèle à 3 niveau.

Modéliser des variances partiellement hétérogènes

Dans notre exemple sur le poids à la naissances des rats, on pouvait penser que les variances ne sont pas hétérogènes entre tous les groupes. Elles pourraient être homogènes pour les groupes ‘Traitements’ qui se distingueraient du groupe ‘Contrôle’. On va donc ajouter une nouvelle variable à la base de données où nous considérons que les groupes “high” et “low” sont identiques et le groupe “control” se différencie.

 rat$trtgrp[rat$Treatment == "Control"] <- "control"
 rat$trtgrp[rat$Treatment == "Low" | rat$Treatment == "High"] <- "traitement"
 str(rat) # constater qu'il y a une nouvelle variable appelée "trtgrp"

## Classes 'tbl_df', 'tbl' and 'data.frame':    322 obs. of  9 variables:
##  $ id       : Factor w/ 322 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ weight   : num  6.6 7.4 7.15 7.24 7.1 6.04 6.98 7.05 6.95 6.29 ...
##  $ sex      : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Litter   : Factor w/ 27 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Lsize    : num  12 12 12 12 12 12 12 12 12 12 ...
##  $ Treatment: Factor w/ 3 levels "Control","High",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "contrasts")= num [1:3, 1:2] 0 1 -1 2 -1 -1
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr  "Control" "High" "Low"
##   .. .. ..$ : NULL
##  $ residus  : num  -0.67 0.884 0.398 0.573 0.301 ...
##   ..- attr(*, "label")= chr "Standardized residuals"
##  $ residus2 : num  -0.884 1.094 0.476 0.699 0.352 ...
##  $ trtgrp   : chr  "control" "control" "control" "control" ...

Note : pour les débutants, il est plus simple de prérarer le fichier de données dans excel par exemple.

Nous pouvons refaire le modèle avec des variances hétérogènes, mais en utilisant la variable trtgrp comme variable permettant d’identité de la variance.

 model.var1 <- lme(weight ~ Treatment + sex + Lsize +Treatment:sex, random = ~ 1 |Litter ,data = rat, method = "REML",
                     weights = varIdent(form = ~1 | Treatment))
 model.var2 <- lme(weight ~ Treatment + sex + Lsize +Treatment:sex, random = ~ 1 |Litter ,data = rat, method = "REML",
                     weights = varIdent(form = ~1 | trtgrp))

A présent, nous pouvons voir si la variance hétérogène sur tous les groupes est un meilleur modèle que les variances hétérogènes uniquement entre les contrôles vs. low/high.

 anova(  model.var1,  model.var2)

##            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## model.var1     1 11 389.0517 430.3300 -183.5259                        
## model.var2     2 10 388.2478 425.7735 -184.1239 1 vs 2 1.196053  0.2741

Comme la probabilité associée à la différence entre les deux modèles est supérieure à 0.05 (i.e., 0.2741), cela signifie que modéliser des variances hétérogènes entre tous les groupes n’améliore pas le modèle. Comme il est plus parcimonieux de modéliser une seule différence (\(Contrôle \neq Traitement\)) plutôt que 3 (i.e., \(Contrôle \neq Haut ;Contrôle \neq Bas ; Bas \neq Haut\))