AVIS IMPORTANT

Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé.

Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le site web de la série d'ateliers R du CSBQ. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés.

Merci de votre compréhension,

Vos coordonnateurs de la série d’ateliers R du CSBQ.

QCBS R Workshops

Cette série de 10 ateliers guide les participants à travers les étapes requises afin de maîtriser le logiciel R pour une grande variété d’analyses statistiques pertinentes en recherche en biologie et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R.

Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil

AVIS IMPORTANT: MISES À JOUR MAJEURES

Mise à jour de mars 2021: Ce wiki n'est plus activement développé ou mis à jour. Le matériel mis à jour pour la série d'ateliers R du CSBQ est maintenant hébergé sur la page GitHub des ateliers R du CSBQ.

Le matériel disponible inclut;

  1. La présentation Rmarkdown pour cet atelier;
  2. Un script R qui suit la présentation.
  3. Le matériel écrit qui accompagne la présentation en format bookdown.

Mise à jour de janvier 2021:

Si vous cherchez l'atelier qui couvre les modèles linéaires généralisés, veuillez vous référer à l'Atelier 6.

Si vous cherchez l'atelier qui couvre les modèles linéaires et généralisés linéaires mixtes, veuillez vous référer à l'Atelier 7.

Atelier 7: Modèles linéaires et généralisés linéaires mixtes

Développé par : Catherine Baltazar, Dalal Hanna, Jacob Ziegler

La section des modèles généralisées linéaires à effets mixtes développé par : Cédric Frenette Dussault, Vincent Fugère, Thomas Lamy, Zofia Taranu

Résumé: Les modèles à effets mixtes permettent aux écologistes de surmonter un certain nombre de limitations liées aux modèles linéaires traditionnels. Dans cet atelier, vous apprendrez à déterminer si vous devez utiliser un modèle à effets mixtes pour analyser vos données. Nous allons vous guider à travers les étapes nécessaires pour utiliser un modèle linéaire mixte, vérifier les suppositions de base et présenter les résultats de votre modèle dans R. Nous allons également nous baser sur ce qui a été appris durant le dernier cours afin d’introduire les modèles linéaires généralisés avec effets mixtes.

Lien vers la nouvelle présentation Rmarkdown

S'il vous plaît essayez-la et dites aux coordonnateurs des ateliers R ce que vous en pensez!

Lien vers l'ancienne présentation Prezi

Télécharger le script de R et les données pour cette leçon:

  1. Données des poissons
  2. Code R
  3. Arabidopsis
  4. Invertébrés
  5. Script GLMM (ce script est pour l'Atelier 6, mais les lignes 446-843 couvrent la section des GLMMs de cet atelier)
  6. glmm_funs (fonctions additionelles pour la section GLMM)

De plus, vous pouvez consulter une présentation plus à jour et plus détaillé sur les modèles mixtes.

Objectifs

  1. Qu’est-ce qu’un modèle linéaire à effets mixtes (MLM) et pourquoi est-ce important?
  2. Comment appliquer des MLMs dans R?
    1. Construction et exploration a priori des modèles et données
    2. Codage et sélection des modèles potentiels
    3. Validation des modèles
    4. Interprétation des résultats et visualisation du modèle
  3. Comment utiliser les MLM quand les suppositions de normailité et d'homogénéité de la variance ne sont pas respectées (GLMMs).

Aperçu

Les données écologiques et biologiques peuvent être complexes et désordonnées. Il existe généralement une structure particulière dans les données (i.e. les observations individuelles ne sont pas toujours indépendantes), les relations entre les variables d'intérêt peuvent différer en fonction de facteurs de regroupement, telles que les espèces, et la plupart du temps de faibles tailles d'échantillons rendent l’ajustement de modèles avec de nombreux paramètres difficile. Les modèles linéaires à effets mixtes (MLM) ont été développés pour aborder ces questions. Ils peuvent être appliqués à un grand nombre de questions écologiques et prennent de nombreuses formes différentes. Dans cet atelier, nous allons utiliser une simple approche interrogative pour apprendre les bases du fonctionnement des MLMs et nous verrons comment les ajuster. Nous terminerons par une réflexion sur l'application des MLMs à un ensemble de questions plus diversifiées.

1. Qu’est-ce qu’un MLM et pourquoi est-ce important?

Les MLMs vous permettent d'utiliser toutes les données que vous avez au lieu d'utiliser des moyennes d'échantillons non indépendants, ils tiennent compte de la structure de vos données (par exemple, quadrats nichés dans les sites eux-mêmes nichés dans les forêts), ils permettent aux relations de varier en fonction de différents facteurs de regroupement (parfois appelés des effets aléatoires) et ils nécessitent moins d’estimation de paramètres que la régression classique, ce qui vous permet d'économiser des degrés de liberté. Mais comment font-ils tout cela? Ces questions vous paraîtront beaucoup plus claires apès avoir lu cette section. Tout d'abord, commençons par se familiariser avec l'ensemble des données.

1.1 Introduction au jeu de données

| Chargez vos données 
# Supprimer commandes antérieures en R
rm(list=ls()) 
 
# Placez tout le matériel de l'atelier dans un même dossier sur votre ordinateur
# Exécutez la ligne de code suivante et utilisez la fenêtre de navigation pour sélectionner QCBS_W6_Data.csv
# le fichier dans le dossier qui contient le matériel de l'atelier
file.choose()
 
# Réglez le répertoire de travail vers le dossier qui contient les fichiers en copiant toute 
# la sortie de la commande file.choose (), à l'exception du nom de fichier R, et collez le tout dans setwd().
 
# Par exemple collez "/Users/ziegljac/Documents/QCBS_R/" -> incluez les guillemets
# NON "/Users/ziegljac/Documents/QCBS_R/Get_Data_Func.R" -> incluez les guillemets 
setwd()
 
# Chargez les bibliothèques et les données utiles
# Si vous n’avez jamais chargé ces bibliothèques avant, vous devrez utiliser la fonction install.packages()
# avant la fonction library()
library(ggplot2)
library(lme4)
library(arm)
library(AICcmodavg)
 
data <- read.csv('QCBS_W6_Data.csv')

Le jeu de données que nous utiliserons porte sur les positions trophiques de poissons. Trois espèces ont été sélectionnées pour l'étude et dix individus par espèce ont été mesurés (longueur corporelle) dans six lacs différents. Voici une représentation visuelle pour vous aider à comprendre tout cela ! Notez : seulement trois individus sont montrés par espèce, mais en réalité il y a 10 individus par espèce.

Une simple question à laquelle vous pourriez répondre avec ce jeu de données est : “est-ce que la position trophique des poissons augmente avec leur taille” ? Au cours de l'atelier, on tentera de répondre à cette question.

DÉFI 1

Pour votre premier défi, vous devez reproduire les graphiques 1 à 3 du code QCBS_W5_LMM.R. Regardez les graphiques et essayez d'obtenir une idée de ce qui se passe. Deux questions clés à se poser sont :

  1. Est-ce qu'on s'attend à ce que, pour toutes les espèces, la position trophique augmente avec la longueur corporelle? Exactement de la même façon?
  2. S'attend-on à ce que ces relations soient pareilles entre les lacs ? Comment pourraient-elles différer ?

Réponse au défi 1

1.2 Comment pourrions-nous analyser ces données?

Nombreuses analyses séparées

Une façon d'analyser ces données est de faire une analyse séparée pour chaque espèce et chaque lac. Voici le graphique de la position trophique en fonction de la taille pour l'espèces 1 dans le lac 1 :

Remarquez que vous devez estimer une ordonnée à l’origine et une pente pour chaque régression (2 paramètres x 3 espèces X 6 lacs = 36 paramètres estimés) et la taille d'échantillon pour chaque analyse est de 10. Il y a peu de chances de détecter un effet à cause de la faible taille d'échantillon et un taux d'erreur augmenté en raison de comparaisons multiples.

Une analyse groupée

Une autre façon d'analyser ces données est de faire une seule analyse en ignorant les variables espèce et lac. Encore une fois, voici le graphique de toutes les données :

Notez que vous avez maintenant une énorme taille d'échantillon et beaucoup moins de paramètres à estimer. Mais qu'en est-il de la pseudoréplication? Les poissons d'un même lac et d'une même espèce sont plus similaires entre-eux que les poissons de lacs et d'espèces différentes. De plus, regardez tout ce bruit dans les données ! Une partie doit être due aux effets de l'espèce et du lac.

Pour notre question, on veut seulement savoir s'il y a un effet général de la longueur corporelle sur la position trophique. On ne cherche pas directement à savoir si cet effet pourrait varier faiblement par espèce en raison de processus biologiques que nous n’avons pas mesurés ou parmi les lacs en raison de variables environnementales non mesurées. Nous voulons simplement prendre en compte cette variation dans le modèle (parfois désignée effets aléatoires).

Construire un MLM

Les MLMs sont un compromis entre séparer et regrouper. Ils:

  1. Estiment une pente et une ordonnée à l'origine pour chaque espèce et chaque lac (séparer), mais en calculant moins de paramètres qu’une régression classique.
  2. Utilisent toutes les données disponibles (regrouper) tout en contrôlant les différences entre les lacs et les espèces (pseudo-réplication).

1.3 Effets fixes et aléatoires

Il y a un débat dans la littérature sur la définition des effets fixes et aléatoires. Il existe plusieurs définitions possibles des effets fixes et aléatoires et nous vous présenterons aujourd'hui celles que nous trouvons plus faciles à appliquer.

Effet fixe

Quand une variable a un effet fixe, les données proviennent de tous les niveaux possibles d'un facteur (variable qualitative). On souhaite émettre des conclusions à propos des niveaux du facteur d'où les données proviennent.

Exemple d'un effet fixe : comparer la concentration de mercure dans les poissons de trois habitats différents. L'habitat est un effet fixe (les trois ont été échantillonnés) et nous sommes intéressés à tirer des conclusions sur les effets de ces trois habitats spécifiques.

Effet aléatoire

Les variables avec un effet aléatoire sont également appelées facteurs aléatoires, car elles sont seulement qualitatives (pas de variable continue). Un effet aléatoire est observé lorsque les données incluent seulement un échantillon aléatoire de tous les niveaux possibles du facteur, qui sont tous d'intérêt. Ils correspondent souvent à des facteurs de regroupement pour lesquels vous souhaitez contrôler l'effet dans votre modèle; vous ne vous intéressez pas à leur effet spécifique sur la variable de réponse.

Exemple d'un effet aléatoire: une étude de la contamination du mercure dans les poissons de lacs de cratères ougandais. Pour des raisons logistiques, vous ne pouvez pas échantillonner tous les lacs de cratères, donc vous échantillonnez seulement huit d'entre eux. Cependant, les poissons d'un lac donné pourrait avoir une sorte de corrélation entre eux (pseudo-corrélation), car ils sont soumis aux mêmes conditions environnementales. Même si vous n'êtes pas intéressé par l'effet de chaque lac spécifiquement, vous devez tenir compte de cette corrélation potentielle avec un facteur aléatoire (lac de cratère) afin de tirer des conclusions sur les lacs de cratères en général.

1.4 Comment fonctionnent les MLMs?

A. Permettre aux intercepts et/ou pentes de varier selon le lac et l'espèce

Permettre aux intercepts et/ou pentes de varier selon certains facteurs (effets aléatoires) signifie simplement que vous supposez qu'ils proviennent d'une distribution normale. La moyenne et l'écart-type de cette distribution sont évalués en fonction de vos données. Les intercepts et pentes les plus probables de cette distribution sont ensuite ajustées par optimisation (ex. maximum de vraisemblance ou maximum de vraisemblance restreint).

Intercepts

Dans le cas d'espèces, seulement la moyenne et l'écart-type de la distribution d’intercepts sont estimés au lieu de trois intercepts pour chaque espèce. La moyenne de cette distribution est le «modèle au niveau de l'espèce». Dans cet exemple, nous n'avons que trois espèces. En général, plus votre facteur comporte de niveaux, plus la moyenne et l'écart-type de la distribution normale seront estimés précisément. Trois niveaux c'est un peu faible, mais plus facile à visualiser! Lorsque vous implémenterez un MLM dans R, notez que l'intercept dans le résumé est l'intercept au niveau des espèces (c.-à-d. la moyenne de tous les intercepts aléatoires).

C'est la même chose pour les lacs : seuls la moyenne et l'écart-type des intercepts des lacs sont estimés au lieu de six intercepts pour chaque lac. Cela économise des degrés de liberté (moins d’estimation de paramètres sont nécessaires).

Pentes

Le même concept s’applique aux pentes qui varient par un facteur donné (effet aléatoire). C’est juste plus difficile à visualiser. Comme dans le cas des espèces, seuls la moyenne et l’écart-type des pentes sont estimés au lieu de trois pentes distinctes. Encore une fois, quand vous implémenterez votre MLM dans R, la pente dans le résumé sera la pente au niveau des espèces.

B. Les intercepts, pentes, et intervalles de confiance associés sont ajustés pour tenir compte de la structure des données

Si une certaine espèce ou un lac est peu représenté (faible échantillon) dans les données (nous avons un design équilibré ici, donc ce n'est pas le cas dans notre exemple), le modèle va accorder plus d'importance au modèle groupé pour estimer l'ordonnée à l'origine et la pente de cette espèce ou de ce lac.

Les intervalles de confiance des intercepts et pentes sont ajustés pour tenir compte de la pseudo-réplication basée sur le coefficient de corrélation intra-classe (CIC) - Combien de variation y a-t-il dans chaque groupe versus entre les groupes ? Ceci détermine votre taille effective de l'échantillon - Une taille d'échantillon ajustée en fonction de la façon dont les données sont corrélées à l’intérieur des groupes.

CIC élevé

Dans ce scénario, le MLM traitera les points provenant du même lac plus comme une moyenne globale, car ils sont fortement corrélés. Par conséquent, la taille effective de l'échantillon sera plus petite, ce qui entraînera des intervalles de confiance plus grands pour la pente et l'intercept.

CIC faible

Dans ce scénario, le MLM traitera les points parvenant du même lac plus indépendamment parce qu’ils sont moins corrélés au sein du groupe que parmi les groupes. Par conséquent, la taille de l'échantillon sera plus grande, ce qui entraînera des intervalles de confiance plus petits pour la pente et l'intercept.

DÉFI 2

Pour votre deuxième défi, répondez aux deux questions suivantes avec vos voisins. Comment le CIC et son intervalle de confiance seront affectés dans ces deux scénarios?

  1. Les positions trophiques des poissons ne varient pas entre les lacs?
  2. Les positions trophiques des poissons sont similaires dans les lacs mais différentes entre les lacs?

Réponse défi 2

2. Comment implémenter un MLM dans R?

Le protocole des modèles mixtes dans R: 

2.1: Construction du modèle a priori et exploration des données
2.2: Coder les modèles potentiels et sélectionner le meilleur modèle
2.3: Valider le modèle
2.4: Interpréter les résultats et visualiser le modèle

2.1 Construction du modèle a priori et exploration des données

Nous voulons déterminer si la position trophique peut être prédite par la longueur corporelle, tout en prenant en compte la variation entre les espèces et les lacs. Donc nous voulons un modèle qui ressemble à ceci:

{PT_ijk} ∼ {Longueur_i} + {Lac_j} + {Espèce_k} + {ε}

où,

{PT_ijk} est la position trophique du poisson (i) du lac (j) et de l’espèce (k)
et
ε sont les résidus du modèle (c. à d. la variation inexpliquée).

Exploration des données

Assurez-vous d'avoir fait le ménage de l'espace de travail (Housekeeping) avant de construire un modèle ! Les données ont-elles la bonne structure ? Utilisez le code suivant pour visualiser la structure et les types de variables au sein de votre jeu de données.

| Exploration des données 
################Section 2#########################
# Exécution d'un modèle mixte en R
# Processus de quatre étapes pour la construction d'un modèle mixte en R
 
# 1) Construction du modèle a priori et exploration des données ####
# i)  Définir un  modèle basé sur une connaissance a priori
# Nous savons que nous voulons construire un modèle qui évalue la relation entre la position trophique et 
# la longueur tout en tenant compte de la variation due au lac et à l'espèce
# Position trophique ~ Longueur + Espèce + Lac
 
# ii)Housekeeping et exploration des données
# Assurez-vous que la structure de vos données soit correcte
str(data)

Sortie

Regardez la distribution des échantillons pour chaque facteur :

| Exploration de données 
# Regardez la distribution des échantillons de chaque facteur pour vérifier
# si le jeu de données est équilibré
table(data$Lake)
table(data$Fish_Species)

Sortie

Ce jeu de donné est parfaitement équilibré, mais les modèles mixtes peuvent analyser les plans expérimentaux non équilibrés (comme c'est souvent le cas en écologie!)

Regardez la distribution des variables continues :

| Exploration de données 
# Regardez la distribution des variables continues
# Transformez si nécessaire (ça évitera des problèmes d’hétérogénéité des résidus du modèle)
hist(data$Fish_Length)
hist(data$Trophic_Pos)

Sortie

Des déviations majeures pourraient causer des problèmes d'hétéroscédasticité. Si c'est nécessaire, appliquez des transformations à vos données. Dans ce cas-ci, les données semblent correctes.

Vérifier la colinéarité entre vos variables explicatives : Vous ne pouvez pas inclure deux variables explicatives colinéaires dans un même modèle, car leurs effets sur la variable réponse seront confondus, c.-à-d. que le modèle ne peut pas indiquer quelle variable colinéaire est responsable de la variation de la variable réponse. Par défaut, le modèle attribuera beaucoup de pouvoir explicatif à la première variable du modèle et peu de pouvoir aux variables qui suivent.

Dans cet exemple, il n’y a pas de risque de colinéarité avec seulement une variable explicative continue. Si vous aviez une autre variable explicative (Var2) et vouliez vérifier la colinéarité, vous pouvez utiliser le code suivant :

| Exploration de données 
# Évaluer la colinéarité entre variables
plot(data)
cor(data$Fish_Length, data$Var2)

DÉFI 3
Quelles mesures supplémentaires aurions-nous pu prendre sur le terrain et qui auraient pu être fortement corrélées avec la longueur corporelle ?

Réponse défi 3

Considérez l'échelle de vos données :
Si deux variables dans le même modèle ont des valeurs se situant sur des échelles très différentes, il est probable que le modèle mixte indique un 'problème de convergence' en essayant de calculer les paramètres. La correction Z standardize les variables et résout ce problème. Elle permet également de mettre toutes vos variables sur la même échelle, même si elles étaient à l'origine de différentes unités :
{z} = ({x} - {mean(x)}) / {sd(x)}

Parce que nos données ont des échelles très différentes, (la longueur est à une échelle plus longue que la position trophique), on applique la correction Z.

| Exploration des données
# Considérez l'échelle de vos données 
# Note : Si deux variables dans le même modèle ont des valeurs se situant sur des échelles très différentes, 
# il est probable que le modèle mixte indique un 'problème de convergence' en essayant de calculer les 
# paramètres. La correction Z standardize les variables et résout ce problème :
# Qu'est-ce qu'une correction de Z ?:  (z = (x - mean(x))/sd(x))
# Longueur corrigée
data$Z_Length<-(data$Fish_Length-mean(data$Fish_Length))/sd(data$Fish_Length)
# Position trophique corrigée
data$Z_TP<-(data$Trophic_Pos-mean(data$Trophic_Pos))/sd(data$Trophic_Pos)

Pour savoir si un modèle mixte est nécessaire pour vos données, vous devez déterminer s'il est important de prendre en compte l'effet aléatoire de facteurs qui pourraient influencer la relation qui vous intéresse (dans notre cas, Lac et Espèce).
Nous pouvons le faire en :

  1. Créant un modèle linéaire sans les facteurs qui pourraient avoir un effet aléatoire
  2. Calculant les résidus de ce modèle linéaire
  3. Produisant un graphique de la valeur des résidus en fonction des niveaux des facteurs potentiellement aléatoires
| Exploration de données 
# Déterminez s'il est important de tenir compte des variations dans les "effets aléatoires" en comparant 
# les résidus d'un modèle linéaire sans les effets aléatoires en fonction des effets aléatoires potentiels
lm.test<-lm(Z_TP~Z_Length, data=data)
lm.test.resid<-rstandard(lm.test)
# Effet de l’espèce
plot(lm.test.resid~ data$Fish_Species, xlab = "Species", ylab="Standardized residuals")
abline(0,0, lty=2)
# Effet du lac
plot(lm.test.resid~ data$Lake, xlab = "Lake", ylab="Standardized residuals")
abline(0,0, lty=2)

Sortie

Pour ce modèle, nous devrions garder les effets aléatoires parce que les résidus standardisés montrent une variation à travers le lac et les espèces.

2.2 Coder les modèles potentiels et sélectionner le meilleur modèle

Notre modèle a priori

{PT_ijk} ∼ {Longueur_i} + {Lac_j} + {Espèce_k} + {ε}

Dans R, on le code ainsi :

La structure de lmer n’est pas intuitive. Les éléments de base de la fonction sont :

REML (Restricted Maximum Likelihood) est la méthode par défaut dans la fonction “lmer”.

Il est à noter que l'estimateur de l’écart-type (sigma) du maximum de vraisemblance (ML) est biaisé d’un facteur (n-2) / n. REML corrige ce biais en faisant un truc (multiplication matricielle de Y telle que la dépendance à X est enlevée). En règle générale , on devrait comparer les modèles d'effets aléatoires nichés avec REML (parce que nous examinons les composantes de la variance et ceux-ci doivent être sans biais), mais lorsque l'on compare des modèles nichés à effets fixes, nous devrions utiliser ML (parce REML fait ce truc avec les matrices X et Y pour corriger le biais).

Mais comment faire si on souhaite aussi que la pente puisse varier?

DÉFI 4
Réécrivez le code suivant de façon à ce que les pentes de la relation de la position trophique en fonction de la longueur corporelle varient par lac et par espèce.

| Coder les MLMs
lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)

Réponse défi 4

Pour déterminer si vous avez construit le meilleur modèle mixte basé sur vos connaissances a priori, vous devez comparer ce modèle a priori à d'autres modèles possibles. Avec le jeu de données sur lequel vous travaillez, il y a plusieurs modèles qui pourraient mieux correspondre à vos données.

DÉFI 5
Faites une liste de 7 modèles alternatifs qui pourraient être construits et comparés à partir de celui-ci: Note - Si nous avions différents effets fixes entre les modèles, nous aurions dû indiquer «REML = FALSE» afin de comparer les modèles avec des méthodes de vraisemblance tels que AIC (voir ci-dessus). Cependant, vous devez rapporter les estimations des paramètres du «meilleur» modèle en utilisant «REML = TRUE».

| Coder les MLMs
lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)

Réponse défi 5

Maintenant que nous avons une liste de modèles potentiels, nous voulons les comparer entre eux pour sélectionner celui (ceux) qui a (ont) le plus grand pouvoir prédictif. Les modèles peuvent être comparés en utilisant la fonction “AICc” provenant du paquet (package) “AICcmodavg”. Le critère d'information d'Akaike (AIC) est une mesure de qualité du modèle pouvant être utilisée pour comparer les modèles. L'AICc corrige pour le biais créé par les faibles tailles d'échantillon quand le AIC est calculé.

Nous allons aussi construire le modèle linéaire de base lm() parce qu'il est toujours utile de voir la variation dans les valeurs de AICc. Pour cette comparaison il est important de changer la méthode à ML (REML=FALSE) parce que lm() n'utilise pas la même méthode d'estimation que lmer(). Par contre, il y a une preuve qui démontre que pour les modèles linéaires de bases les résultats de la méthode des moindres carrés (Least squares) est équivalente au résultats de la méthode ML.

| Comparer les MLMs
# 2) Coder les modèles potentiels et sélectionner le meilleur modèle ####
# i) Coder les modèles potentiels
# Liste de tous les modèles potentiels -->
# Note: vous pouvez choisir de ne pas coder ceux qui n'ont pas de sens biologique.
# Construisez aussi le modèle lm() pour voir la variation dans les valeurs de AICc, mais changez la 
# méthode à ML (REML=FALSE) parce que lm() n'utilise pas la même méthode d'estimation que lmer().
# Modèle linéaire sans effets aléatoires
M0<-lm(Z_TP~Z_Length,data=data)
# Modèle complet avec différents intercepts
M1<-lmer(Z_TP~Z_Length + (1|Fish_Species) + (1|Lake), data=data, REML=FALSE)
# Modèle complet avec différents intercepts et pentes
M2<-lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species) + (1+Z_Length|Lake), data=data, REML=FALSE)
# Aucun effet Lac, intercept aléatoire seulement
M3<-lmer(Z_TP~Z_Length + (1|Fish_Species), data=data, REML=FALSE)
# Aucun effet Espèce, intercept aléatoire seulement
M4<-lmer(Z_TP~Z_Length + (1|Lake), data=data, REML=FALSE)
# Aucun effet Lac, intercept et pente aléatoires
M5<-lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species), data=data, REML=FALSE)
# Aucun effet Espèce, intercept et pente aléatoires
M6<-lmer(Z_TP~Z_Length + (1+Z_Length|Lake), data=data, REML=FALSE)
# Modèle complet avec intercepts et pentes qui variant par Lac
M7<-lmer(Z_TP~Z_Length + (1|Fish_Species) + (1+Z_Length|Lake), data=data, REML=FALSE)
# Modèle complet avec intercepts et pentes qui variant par Espèce
M8<-lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species) + (1|Lake), data=data, REML=FALSE)
 
 
# ii) Comparer les modèles en comparant les valeurs AICc 
# Calculer les valeurs AICc pour chaque modèle
AICc<-c(AICc(M0), AICc(M1), AICc(M2), AICc(M3), AICc(M4), AICc(M5), AICc(M6), AICc(M7), AICc(M8))
# Mettre des valeurs dans une table pour faciliter la comparaison
Model<-c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8")
AICtable<-data.frame(Model=Model, AICc=AICc)
AICtable
# M8 a la plus faible valeur AICc donc le meilleur modèle
# M2 est également un bon modèle, mais tous les autres modèles ne sont pas aussi bons.

Sortie

Le modèle avec un AICc plus bas a le plus grand pouvoir prédictif considérant les données. Certains disent que si deux modèles sont à plus ou moins 2 unitées d'AICc de différence, leurs pouvoirs prédictifs sont équivalents. Dans notre cas, on peut regarder de plus près M8 et M2, mais tous les autres ont des AICc tellement plus élevés qu'on peut exclure la possibilité qu'ils soient les meilleurs modèles pour nos données.

| Recoder les MLMs
# Une fois que les meilleurs modèles sont sélectionnés il faut remettre REML=TRUE
M8<-lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species) + (1|Lake), data=data, REML=TRUE)
M2<-lmer(Z_TP~Z_Length + (1+Z_Length|Fish_Species) + (1+Z_Length|Lake), data=data, REML=TRUE)

DÉFI 6
Prenez 2 minutes avec votre voisin pour étudier la structure du modèle M2. Comment diffère-t-elle de M8 d'un point de vue biologique? Pourquoi n'est-il pas surprenant que sa valeur de AICc soit la deuxième meilleure?

Réponse défi 6

2.3 Validation du modèle

Pour vérifier la supposition d'homogénéité, il faut faire un graphique des valeurs prédites en fonction des valeurs résiduelles.

| Validation du modèle 
#3) Vérification des suppositions du modèle ####
# Vérification pour M8
# A. Vérifiez l'homogénéité : graphique des valeurs prédites vs valeurs résiduelles
E1 <- resid(M8)
F1 <- fitted(M8)
plot(x = F1, 
     y = E1, 
     xlab = "Fitted Values",
     ylab = "Normalized residuals")
abline(h = 0, lty = 2)

Sortie

L'étendue similaire des résidus suggère que le modèle est adéquat pour bien modéliser nos données.

Pour vérifier la supposition d'indépendance, il faut faire un graphique des résidus en fonction de chaque covariable du modèle :

| Validation du modèle 
# B. Vérifiez l’indépendance :
# i. graphique des résidus VS chaque covariable du modèle 
# Longueur corporelle des poissons
plot(x = data$Z_Length, 
     y = E1, 
     xlab = "Z Length",
     ylab = "Normalized residuals")
abline(h = 0, lty = 2)
# Note: Les regroupements de données sont dus à la structure des données, où des poissons de seulement 5 
# classes de taille  (grandes, petites, et trois groupes entre les deux) étaient capturés.
 
# Espèce
boxplot(E1 ~ Fish_Species,   
        ylab = "Normalized residuals",
        data = data, xlab = "Species")
abline(h = 0, lty = 2)
 
# Lac
boxplot(E1 ~ Lake,   
        ylab = "Normalized residuals",
        data = data, xlab = "Lake")
abline(h = 0, lty = 2)

Sortie

L'étendue similaire au dessus et sous zéro indique qu'il n'y a pas de problème d'indépendance avec cette variable.

Idéalement, vous devriez aussi faire l'analyse ci-dessus pour chaque covariable non inclus dans votre modèle. Si vous observez des patrons dans ces graphiques, vous saurez qu'il y a de la variation dans votre jeu de données qui pourrait être expliquée par ces covariables et vous devriez considérer d'inclure ces variables dans votre modèle. Puisque dans notre cas, nous avons inclus toutes les variables mesurées dans notre modèle, nous ne pouvons pas faire cette étape.

Il est également important de vérifier la normalité des résidus. Des résidus suivant une distribution normale indiquent que le modèle n'est pas biaisé.

| Validation du modèle 
#D. Vérifier la normalité : histogramme
hist(E1)

Sortie

2.4 Interpréter les résultats et les visualiser graphiquement

Vous pouvez voir le résumé du modèle à l'aide de :

| Interprétation des résultats
# Vérifiez le résumé du modèle
# Cela vous permet d'avoir une idée de la variance expliquée
# par les différentes composantes du modèle et la «significativité» des effets fixes 
summary(M8)

Sortie

La sortie est divisée en descriptions des effets aléatoires (ce qui peut varier en fonction de la distribution normale) et les effets fixes (ce que nous estimons comme pour une régression classique) :

Pour déterminer si la pente, et donc l'effet de la longueur sur la position trophique, est significativement différente de zéro, vous devez d'abord calculer l'intervalle de confiance (IC) du paramètre de la pente (estimation pour Z_Length dans la section des effets fixes = 0,4223). CI = l'erreur-type de l'estimation x 1,96 plus ou moins l'estimation du paramètre. Si l’IC inclut zéro, la pente n’est pas significativement différente de zéro au seuil de 0,05.

DÉFI 7
a) Quelle est la pente et l'intervalle de confiance de la variable Z_Length du le modèle M8?

b) Est-ce que la pente de Z_Length est significativement différente de 0?

Réponse défi 7

Pour mieux visualiser les résultats d'un modèle mixte, les différentes ordonnées à l'origine et pentes générées par le modèle peuvent être représentées dans des figures. Nos coefficients du modèle au niveau du groupe (aka: “coefs”, dans ce cas seulement un intercept et une pente) se trouvent dans le résumé du modèle dans la section des effets fixes. Les “coefs” pour chacun des niveaux du modèle (dans ce cas: Lac et Espèces) qui ont été ajustés à une distribution normale peuvent être obtenus en utilisant la fonction “coef ()”.

Deux façons de visualiser ces données sont :

  1. Montrer le modèle au niveau du groupe (toutes les données groupées)
  2. Montrer le modèle au niveau de l’espèce ou du lac

1. Pour montrer le modèle au niveau du groupe :

Obtenir les paramètres d'intérêts

et tracer les données avec le modèle superposée

| Visualiser le modèle 
# Visualiser les résultats du modèle ####
# Il existe plusieurs façons de visualiser les résultats d'un modèle mixte, qui font tous appel au coefficients générés par le modèle.
# La première étape est d'obtenir les coefficients du modèle afin de les ajouter aux figures
coef(M8)
# Maintenant, mettez les coefs dans un tableau pour les rendre plus faciles à manipuler
Lake.coef <- as.data.frame(coef(M8)$Lake)
colnames(Lake.coef) <- c("Intercept","Slope")
Species.coef <- as.data.frame(coef(M8)$Fish_Species)
colnames(Species.coef) <- c("Intercept","Slope")
 
# Graphique 1 – toutes les données
# Graphique qui inclut toutes les données
plot <- ggplot(aes(Z_Length,Z_TP),data=data)
Plot_AllData <- plot + geom_point() + xlab("Length (mm)") + ylab("Trophic Position") + labs(title="All Data") + fig
# Ajoutez un abline avec l'intercept et la pente de la relation entre la longueur et la position trophique
# Notez que vous pouvez obtenir l’origine et la pente du facteur fixe directement à partir du résumé du modèle
summary(M8)
Plot_AllData + geom_abline(intercept = -0.0009059, slope =0.4222697)

Sortie

2. Montrer les modèle au niveau de l’espèce ou du lac:

Obtenir les paramètres d'intérêts

et tracer les données avec le modèle superposé

| Visualiser le modèle
# Graphique 2 - Par Espèce 
# Colorez les données par espèce
Plot_BySpecies<-plot + geom_point(aes(colour = factor(Fish_Species)), size = 4) + xlab("Length (mm)") + ylab("Trophic Position") + labs(title="By Species") + fig
# Ajoutez les lignes de régression pour chaque espèce
Plot_BySpecies + geom_abline(intercept = Species.coef[1,1], slope =Species.coef[1,2], colour="coral2") + geom_abline(intercept = Species.coef[2,1], slope =Species.coef[2,2], colour = "green4") + geom_abline(intercept = Species.coef[3,1], slope =Species.coef[3,2], colour="blue1")
 
# Graphique 3 – Par Lac 
# Colorez les données par lac
Plot_ByLake<-plot + geom_point(aes(colour = factor(Lake)), size = 4) + xlab("Length (mm)") + ylab("Trophic Position") + labs(title="By Lake") + fig
# Ajouter les lignes de régression avec les intercepts spécifiques à chaque lac
Plot_ByLake + geom_abline(intercept = Lake.coef[1,1], slope =Lake.coef[1,2], colour="coral2") + geom_abline(intercept = Lake.coef[2,1], slope =Lake.coef[2,2], colour="khaki4") + geom_abline(intercept = Lake.coef[3,1], slope =Lake.coef[3,2], colour="green4") + geom_abline(intercept = Lake.coef[4,1], slope =Lake.coef[4,2], colour="darkgoldenrod") + geom_abline(intercept = Lake.coef[5,1], slope =Lake.coef[5,2], colour="royalblue1") + geom_abline(intercept = Lake.coef[6,1], slope =Lake.coef[6,2], colour="magenta3")

Sortie

Exercice de réflexion

Les modèles mixtes sont très utiles pour prendre en compte la structure complexe des données en écologie tout en permettant de ne pas perdre beaucoup de degrés de liberté.

Nous avons couvert seulement une petite partie de ce que les MLMs peuvent faire. Ci-dessous vous trouverez quelques autres exercices avec des structures de données semblables aux données de l'atelier et deux livres qui détaillent l'utilité des MLMs.

DÉFI 8
Situation : Vous avez récolté des estimés de biodiversité dans 1000 quadrats qui sont dans 10 différents sites et qui sont également dans 10 forêts différentes (i.e. 10 quadrats par site et 10 sites par forêt). Vous avez de plus mesuré la productivité dans chaque quadrat. Vous cherchez à savoir si la productivité est un bon prédicteur de la biodiversité.

Quel modèle mixte pourriez-vous utiliser pour ce jeu de données?

Réponse défi 8

DÉFI 9
Situation : Vous avez récolté 200 poissons dans 12 différents sites distribués également dans 4 habitats différents qui se retrouvent dans un même lac. Vous avez mesuré la longueur de chaque poisson et la quantité de mercure dans ses tissus. Vous cherchez surtout à savoir si l'habitat est un bon prédicteur de la concentration en mercure.

Quel modèle mixte pourriez-vous utiliser pour ce jeu de données?

Réponse défi 9

Les modèles linéaires généralisés mixtes (GLMM)

On vient de voir que les variables réponses binaires et d'abondance peuvent être modélisées de façon plus appropriée en laissant tomber la supposition que les résidus doivent suivre une distribution normale en spécifiant une distribution pour les résidus (par exemple Poisson ou binomiale négative). Toutefois, que devrions-nous faire s'il existe une structure nichée dans les données ou des corrélations entre les observations qui causent certaines observations à être plus semblables les unes aux autres que les observations échantillonnées dans différents sites ou points dans le temps ? Comme nous l'avons vu dans l'atelier 6, cette non-indépendance peut être expliquée en ajoutant des termes à effets aléatoires dans le modèle. Rappelez-vous que les intercepts et / ou pentes peuvent varier en fonction des effets aléatoires (facteurs de regroupement) et que l’écart entre ces pentes et / ou intercepts suit une distribution normale. Les effets aléatoires sont également nommés “estimations de rétrécissement” parce qu'ils représentent une moyenne pondérée des données et de l’effet global (effet fixe). La quantité de rétrécissement vers l'effet global est plus sévère si la variabilité intra-groupe est large par rapport à la variabilité inter-groupe (i.e. quand nous avons moins confiance en l'estimation aléatoire, en raison de la faible taille de l'échantillon ou d'une haute variabilité intra-groupe, l'estimation est 'tirée' vers l'effet fixe).

Les propriétés des modèles linéaires mixtes (atelier 6) qui incorporent des effets aléatoires et celles des modèles linéaires généralisés (vu ci-dessus), qui permettent de spécifier une distribution statistique autre que la distribution normale, peuvent être combinées dans le cadre des modèles linéaires mixtes généralisés (GLMMs). L'extension des GLM qui spécifie cette structure supplémentaire suit des étapes similaires à celles présentées dans l'atelier sur les modèles linéaires mixte.

Pour donner un bref aperçu des GLMM, nous allons voir une étude de cas présentée par Bolker et al. (2009) et par la suite Bolker et al. (2011) où les observations échantillonnées dans différentes populations ont créé une structure (ou manque d'indépendance) dans le jeu de données. Plus précisément, les auteurs ont évalué comment l'abondance d'Arabidopsis thaliana (arabette des dames) varie en fonction de variables environnementales (disponibilité d'éléments nutritifs et herbivorie) et en fonction de leur population d'origine et / ou de leur génotype.

En utilisant l’ensemble des données d'Arabidopsis, nous allons maintenant examiner l'effet des niveaux de nutriments, d'herbivorie et leur interaction (effets fixes) sur la production de fruits d'Arabidopsis thaliana et la variabilité de ces relations à travers différentes populations et / ou génotypes (effets aléatoires).

| Charger les données 
# Chargez et affichez le jeu de données
dat.tf <- read.csv("Banta_TotalFruits.csv")
str(dat.tf)
 
# X         numéro d'observation
# reg       facteur pour l'une des trois régions; Pays-Bas, l'Espagne ou la Suède
# popu      facteur avec un niveau pour chaque population (effet aléatoire)
# gen       facteur avec un niveau pour chaque génotype (effet aléatoire)
# rack      facteur de nuisance pour l'un des deux racks placés dans la serre
# nutrient  facteur à deux niveaux précisant une faible (valeur = 1) ou haute (valeur = 8) concentration 
#           de nutriments (effet fixe)
# amd       facteur à deux niveaux précisant l'absence ou la présence d'herbivorie (lésions sur le méristème apical)
#           (effet fixe)
# status    facteur de nuisance pour la méthode de germination
# total.fruits  variable réponse; nombre entier indiquant le nombre de fruits par plante
 
# 2-3 génotypes imbriqués dans chacune des neuf populations
table(dat.tf$popu,dat.tf$gen)
 
# Entretien : changer nombres entiers en facteurs, rajuster les niveaux d'herbivorie (amd) et renommer les 
#             niveaux de nutriments
dat.tf <- transform(dat.tf,
                    X=factor(X),
                    gen=factor(gen),
                    rack=factor(rack),
                    amd=factor(amd,levels=c("unclipped","clipped")),
                    nutrient=factor(nutrient,label=c("Low","High")))
 
# Installer / charger les librairies
if(!require(lme4)){install.packages("lme4")}
require(lme4)
if(!require(coefplot2)){install.packages("coefplot2",repos="http://www.math.mcmaster.ca/bolker/R",type="source")}
require(coefplot2)     
if(!require(reshape)){install.packages("reshape")}
require(reshape) 
if(!require(ggplot2)){install.packages("ggplot2")}
require(ggplot2)
if(!require(plyr)){install.packages("plyr")}
require(plyr)
if(!require(gridExtra)){install.packages("gridExtra")}
require(gridExtra)
if(!require(emdbook)){install.packages("emdbook")}
require(emdbook)
source("glmm_funs.R")

La structure du jeu de données

Lorsque nous examinons l'interaction entre les éléments nutritifs et l'herbivorie (clipping) à travers les neuf populations différentes, nous constatons que le nombre de fruits est toujours plus élevé dans le traitement élevé (High) en nutriments. L'effet de clipping est plus faible, même que dans certaines populations nous notons une pente négative.

facet par popu
ggplot(dat.tf,aes(x=amd,y=log(total.fruits+1),colour=nutrient)) +
  geom_point() + 
  stat_summary(aes(x= as.numeric(amd)),fun.y=mean,geom="line") +
  theme_bw() + theme(panel.margin=unit(0,"lines")) +
  scale_color_manual(values=c("#3B9AB2","#F21A00")) + # de la palette Wes Anderson Zissou
  facet_wrap(~popu)

Un graphique similaire peut être fait pour les 24 génotypes différents (en changeant facet_wrap(~gen)).

Choisir une distribution pour les résidus

La variable réponse représente des données d'occurrences, donc nous devons choisir une distribution de Poisson pour modéliser cette variable. Rappelez-vous qu’une propriété importante de la distribution de Poisson est que la variance est égale à la moyenne. Cependant, comme nous le verrons ci-dessous, la variance de chaque groupe augmente beaucoup plus rapidement que la moyenne attendue sous la distribution de Poisson.

| Hétérogénéité entre les groupes 
# Créez de nouvelles variables qui représentent toutes les combinaisons de 
# nutriments x facteur clipping x facteur aléatoire
dat.tf <- within(dat.tf,
{
  # génotype x nutriments x clipping
  gna <- interaction(gen,nutrient,amd)
  gna <- reorder(gna, total.fruits, mean)
  # population x nutriments x clipping
  pna <- interaction(popu,nutrient,amd)
  pna <- reorder(pna, total.fruits, mean)
})
 
 
# Boîte à moustaches du total de fruits vs. nouvelle variable (génotype x nutriments x clipping)
ggplot(data = dat.tf, aes(factor(x = gna),y = log(total.fruits + 1))) +
   geom_boxplot(colour = "skyblue2", outlier.shape = 21, outlier.colour = "skyblue2") + 
   theme_bw() + theme(axis.text.x=element_text(angle=90)) + 
   stat_summary(fun.y=mean, geom="point", colour = "red") 
 
# Boîte à moustaches du total des fruits vs. nouvelle variable (population x nutriments x clipping)
ggplot(data = dat.tf, aes(factor(x = pna),y = log(total.fruits + 1))) +
  geom_boxplot(colour = "skyblue2", outlier.shape = 21, outlier.colour = "skyblue2") + 
  theme_bw() + theme(axis.text.x=element_text(angle=90)) + 
  stat_summary(fun.y=mean, geom="point", colour = "red")

Tel qu'illustré ci-dessus, il existe une importante hétérogénéité parmi la variance de chaque groupe même lorsque la variable réponse est transformée. Nous notons également que certains groupes ont une moyenne et variance de zéro.

Pour identifier la famille de distribution la plus appropriée , nous pouvons examiner un graphique de diagnostic de la variance de chaque groupe par rapport à leurs moyennes. Nous présentons un exemple ci-dessous pour le regroupement par génotype x nutriments x herbivore (clipping).

  1. Si nous observons une relation linéaire entre la variance et la moyenne avec une pente = 1, une famille de Poisson serait appropriée,
  2. Si nous observons une relation moyenne-variance linéaire avec une pente> 1 (c. Var = φµ où φ > 1), la famille quasi-Poisson (tel que présenté ci-dessus) doit être utilisée,
  3. Enfin, une relation quadratique entre la variance et la moyenne (c. Var = µ(1 + α

) ou µ(1 + µ/k)), est caractéristique des données surdispersées résultant d'une hétérogénéité sous-jacente entre les échantillons. Dans ce cas, la distribution binomiale négative (Poisson-gamma) serait plus appropriée.

En examinant la figure ci-dessus, nous constatons qu'une distribution quasi-Poisson linéaire peut être meilleure que la distribution binomiale négative, mais nous aurions besoin de modélisation supplémentaire pour le confirmer.

Un GLMM avec une distribution de Poisson

On crée un GLMM en utilisant la fonction glmer() de la librairie lme4. On spécifie le modèle avec un intercepte aléatoire pour les facteurs génotype et population. Nous incluons les variables de nuisance (rack et la méthode de germination) comme effets fixes. Compte tenu de la relation entre la moyenne et la variance que nous avons observée ci-dessus, nous aurons probablement besoin d'un modèle avec surdispersion, mais nous allons commencer par un modèle avec une distribution de Poisson.

| Poisson GLMM
mp1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
               (1|popu)+
               (1|gen),
             data=dat.tf, family="poisson")

Nous pouvons vérifier la surdispersion en utilisant la fonction overdisp_fun fournie par Bolker et al. (2011) qui divise les résidus de Pearson par les degrés de liberté des résidus et vérifie si les valeurs diffèrent de façon significative (i.e. si le rapport entre la déviance et les degrés de liberté est significativement différent de 1). Une faible valeur de p indique que les données sont surdispersées (i.e. le ratio est significativement différent de 1).

| Vérification de la surdispersion 
# Surdispersion?
overdisp_fun(mp1)
 
# On peut aussi l’estimer en divisant la déviance résiduelle par les degrés de liberté des 
# résidus
summary(mp1)
# déviance résiduelle = 18253.7 and dl resid = 616

Le ratio est » 1, donc comme on s'y attendait, nous devons utiliser une distribution différente où la variance augmente plus rapidement que la moyenne, tels que la distribution binomiale négative.

Un GLMM avec un distribution binomiale négative

Rappelez-vous que la distribution binomiale négative (ou Poisson-gamma) satisfait la supposition que la variance est proportionnelle au carré de la moyenne .

| NB
# Note : Ce modèle converge si vous utilisez la version 3.0.2 de R, mais peut ne pas converger avec des 
# versions plus récentes. Si vous avez des problèmes de convergence, essayez le code suivant avec la
# version 3.0.2.
 
mnb1 <- glmer.nb(total.fruits ~ nutrient*amd + rack + status + 
               (1|popu)+
               (1|gen),
             data=dat.tf, control=glmerControl(optimizer="bobyqa"))
# Le nouvel argument «control», bien qu’au-delà de la portée de cet atelier, spécifie la façon dont nous
# optimisons les valeurs des paramètres (c. en prenant la dérivée d'une fonction ou en procédant par itération).
# Si ce n’est pas possible de prendre la dérivée de la fonction, un algorithme itératif comme bobyqa 
# (Bound Optimization BY Quadratic Approximation) est utilisé.
 
# Surdispersion?
overdisp_fun(mnb1)

Le rapport est maintenant beaucoup plus près de 1 (bien que la valeur p est encore inférieur à 0,05).

Un GLMM avec une distribution "Poisson-lognormal"

Une autre option que nous n’avons pas encore décrit est la distribution Poisson-lognormal. Cette approche aborde le problème de la surdispersion en appliquant un effet aléatoire au niveau de l'observation (voir Harrison (2014)), et modèle la variation Poisson supplémentaire de la variable de réponse en utilisant un effet aléatoire avec un niveau de chaque observation unique.

Un modèle Poisson-lognormal peut être construit en modélisant les εi avec une distribution a priori lognormal. Une distribution Poisson-lognormal avec une moyenne µ et une variance a priori lognormal σ2 est donné par :

var(y) = µ + µ2 [exp(σ2) - 1]

En revanche, nous avons vu que la distribution binomiale négative (Poisson-gamma) était donné par:

var(y) = µ + µ2/k

En générale, la variance σ2 dans la distribution Poisson-lognormal dépendra du niveau de regroupement que nous sélectionnons (i.e. au niveau individuel, génotype ou de la population). Donc, le modèle de Poisson-lognormal nous permet d'assigner l'agrégation observée à différentes sources d'hétérogénéité. Ici, afin d'évaluer un effet aléatoire au niveau de l'observation, nous allons l'évaluer au niveau individuel.

| PL GLMM
mpl1 <- glmer(total.fruits ~ nutrient*amd + rack + status + 
               (1|X) +
               (1|popu)+
               (1|gen),
             data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))
 
overdisp_fun(mpl1)

Nous avons réussi! Le rapport entre la déviance et les degrés de liberté est maintenant conforme avec notre critère (en fait, il est plus petit 1).

Intercepts aléatoires

Maintenant que nous avons la distribution d'erreur appropriée, nous pouvons tester l'importance des interceptes aléatoires (pour population et génotype) en comparant des modèles nichés avec et sans les effets aléatoires d'intérêt en utilisant soit:

  1. l'approche théorique d'information (tel que le Critère d'Information d'Akaike; AIC), qui, comme nous l'avons vu dans l'atelier 6 examine plusieurs hypothèses concurrentes (modèles) simultanément pour identifier le modèle avec le pouvoir prédictif le plus élevé compte tenu des données. Nous allons encore une fois utiliser le AICc pour corriger les petites tailles d'échantillon.
  2. l'approche fréquentiste (test d'hypothèse nulle traditionnelle), où deux modèles nichés sont comparés en utilisant le tests de rapport de vraisemblance de la fonction anova(). Il est important de noter qu’avec cette approche, nous testons une hypothèse nulle d'une variance de zéro, mais étant donné que nous ne pouvons pas avoir un écart négatif, nous testons le paramètre à la limite de sa région réalisable. Par conséquent, la valeur de p rapporté est environ le double de ce qu'elle devrait être (c. nous avons tronquée la moitié des valeurs possibles ; celles qui tombent en dessous de 0).
Évaluation des intercepts aléatoires
summary(mpl1)$varcor
 
# popu seulement
mpl1.popu <- glmer(total.fruits ~ nutrient*amd + rack + status + 
                     (1|X) +
                     (1|popu), 
                     data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))
 
# gen seulement
mpl1.gen <-glmer(total.fruits ~ nutrient*amd + rack + status + 
                   (1|X) +
                   (1|gen), 
                   data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))
 
# Approche AICc
ICtab(mpl1, mpl1.popu, mpl1.gen, type = c("AICc"))
#            dAICc df
# mpl1       0.0   10
# mpl1.popu  2.0   9 
# mpl1.gen   16.1  9 
 
# Approche fréquentiste (Likelihood Ratio Test)
anova(mpl1,mpl1.popu)
# Data: dat.tf
# Df         AIC  BIC      logLik  deviance  Chisq   Chi    Df    Pr(>Chisq)  
# mpl1.popu  9    5017.4   5057.4  -2499.7   4999.4                           
# mpl1       10   5015.4   5059.8  -2497.7   4995.4  4.0639  1    0.04381 *
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
anova(mpl1,mpl1.gen)
# Df        AIC  BIC     logLik deviance  Chisq    Chi     Df   Pr(>Chisq)    
# mpl1.gen  9    5031.5  5071.5 -2506.8   5013.5                             
# mpl1      10   5015.4  5059.8 -2497.7   4995.4   18.177  1   2.014e-05 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Notez que le modèle sans l'intercept aléatoire pour génotype est à moins de deux unités AICc du modèle complet, ce qui indique que les deux modèles sont tout aussi plausibles (i.e. peu de soutien pour inclure un intercept aléatoire pour le génotype). Toutefois, lorsque nous utilisons la comparaison de vraisemblance de modèles nichés (anova()), et corrigeons pour les valeurs p gonflés par un facteur de 2, nous trouvons que dans les deux cas, p <0,05 et donc le modèle avec les deux effets aléatoires (mpl1) est sélectionné.

Représentation graphique des paramètres du modèle

Une représentation graphique des paramètres du modèle peut être obtenue en utilisant la fonction coefplot2. Par exemple, pour afficher les paramètres de variance de nos trois intercepts aléatoires:

Coefplot2
# Paramètres de la variance
coefplot2(mpl1,ptype="vcov",intercept=TRUE,main="Random effect variance")
 
# Effets fixes
coefplot2(mpl1,intercept=TRUE,main="Fixed effect coefficient")

Notez que les barres d'erreur ne sont visibles que pour les effets fixes parce glmer ne nous donne pas d'informations sur l'incertitude des paramètres de variance.

Graphiques des intercepts aléatoires

Nous pouvons également extraire les prédictions des effets aléatoires en utilisant la fonction ranef et les tracer en utilisant la fonction dotplot, qui crée un graphique à deux facettes pour chaque effet aléatoire (popu et gen). Notez : la fonction grid.arrange est utilisée pour supprimer l'effet aléatoire au niveau de l'observation (i.e. (1 | X)).

dotplot
pp <- list(layout.widths=list(left.padding=0, right.padding=0),
           layout.heights=list(top.padding=0, bottom.padding=0))
r2 <- ranef(mpl1,condVar=TRUE)
d2 <- dotplot(r2, par.settings=pp)
grid.arrange(d2$gen,d2$popu,nrow=1)

Le graphique indique un peu de variabilité régionale parmi les populations où les populations espagnoles (SP) ont des valeurs plus élevées que les populations suédoises (SW) et néerlandais (NL). La différence entre les génotypes semble largement induite par le génotype 34.

Pentes aléatoires

Section bonus: Pentes aléatoires

Modèle final

Nous terminons cette section sur les GLMMs avec l'évaluation des effets fixes. Nous décrivons d'abord les modèles potentiels et les comparons en utilisant AICc.

Évaluation des effets fixes avec AICc 
mpl2 <- update(mpl1, . ~ . - rack) # modèle sans rack
mpl3 <- update(mpl1, . ~ . - status) # modèle sans status
mpl4 <- update(mpl1, . ~ . - amd:nutrient) # modèle sans l’ interaction amd:nutrient 
ICtab(mpl1,mpl2,mpl3,mpl4, type = c("AICc"))

Nous pouvons aussi utiliser les fonctions drop1 et dfun, où dfun convertit les valeurs AIC retournées par drop1 en valeurs ΔAIC (produisant un tableau comparable à ICtab ci-dessus).

Évaluation des effets fixes avec drop1
(dd_LRT <- drop1(mpl1,test="Chisq"))
# Model:
#   total.fruits ~ nutrient * amd + rack + status + (1 | X) + (1 | popu) + (1 | gen)
#              Df    AIC   LRT   Pr(Chi)    
# <none>          5015.4                     
# rack          1 5070.5 57.083 4.179e-14 ***
# status        2 5017.0  5.612   0.06044 .  
# nutrient:amd  1 5016.8  3.444   0.06349 .  
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
(dd_AIC <- dfun(drop1(mpl1)))
#               Df   dAIC
# <none>           0.000
# rack          1 55.083
# status        2  1.612
# nutrient:amd  1  1.444

Il y a un fort effet de rack (un changement de AIC de 55 si on enlève cette variable), alors que les effets de « status » et de l'interaction sont faibles (différence AIC moins de 2). Nous pouvons donc commencer par enlever l'interaction non significative afin de tester les effets principaux de nutriments et d’herviborie (clipping).

Supprimer l’interation
mpl2 <- update(mpl1, . ~ . - and:nutrient)
 
# Avec AICc:
mpl3 <- update(mpl2, . ~ . - rack) # modèle sans rack ou l’interaction
mpl4 <- update(mpl2, . ~ . - status) # modèle sans status ou l’interaction
mpl5 <- update(mpl2, . ~ . - nutrient) # modèle sans nutrient ou l’interaction
mpl6 <- update(mpl2, . ~ . - amd) # modèle sans clipping ou l’interaction
ICtab(mpl2,mpl3,mpl4,mpl5,mpl6, type = c("AICc"))
#        dAICc  df
# mpl2   0.0    9 
# mpl4   1.2    7 
# mpl6   10.2   8 
# mpl3   54.2   8 
# mpl5   135.6  8 
 
# Ou, avec drop1:
(dd_LRT2 <- drop1(mpl2,test="Chisq"))
# Model:
#   total.fruits ~ nutrient + amd + rack + status + (1 | X) + (1 | popu) + (1 | gen)
#            Df    AIC     LRT   Pr(Chi)    
#   <none>      5016.8                      
#   nutrient  1 5152.5 137.688 < 2.2e-16 ***
#   amd       1 5027.0  12.218 0.0004734 ***
#   rack      1 5071.0  56.231 6.443e-14 ***
#   status    2 5018.1   5.286 0.0711639 .  
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
(dd_AIC2 <- dfun(drop1(mpl2)))
#           Df    dAIC
# <none>        0.000
# nutrient  1 135.688
# amd       1  10.218
# rack      1  54.231
# status    2   1.286
 
summary(mpl2)

Les effets de nutriments et d’herbivorie sont forts (grand changement d'AIC de 135,6 (mpl5) et 10,2 (mpl6) si nutriments ou clipping sont supprimés respectivement). Notre modèle final inclut l’effet fixe des nutriments (fortement positif) et l’effet fixe d’herbivorie (négatifs), la variable nuisance de rack, l’effet aléatoire au niveau de l'observation (1 | X) pour tenir compte de la surdispersion et la variation de fruits par population et génotype .

DÉFI 10

En utilisant l’ensemble de données inverts (temps de développement larvaire (PLD) de 74 espèces d'invertébrés et vertébrés marins élevés à différentes températures et temps), répondez aux questions suivantes:

  1. Quel est l'effet du type d'alimentation et du climat (effets fixes) sur PLD (variable réponse)?
  2. Est-ce que cette relation varie selon les taxons (effets aléatoires)?
  3. Quelle est la meilleure famille de distribution pour ces données?
  4. Une fois que vous avez déterminé la meilleure famille de distribution, réévaluez vos effets aléatoires et fixes.

Défi 10: Solution 1

Défi 10: Solution 2

Défi 10: Solution 3

Défi 10: Solution 4

Ressources additionnelles

Livres de référence pour les MLMs

Gelman, A., and Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models (Cambridge University Press).

Zuur, A., Ieno, E.N., Walker, N., Saveliev, A.A., and Smith, G.M. (2009). Mixed effects models and extensions in ecology with R (Springer).

Ressources pour les GLMMs

Livres

B. Bolker (2009) Ecological Models and Data in R. Princeton University Press.
A. Zuur et al. (2009) Mixed Effects Models and Extensions in Ecology with R. Springer.

Articles

Harrison, X. A., L. Donaldson, M. E. Correa-Cano, J. Evans, D. N. Fisher, C. E. D. Goodwin, B. S. Robinson, D. J. Hodgson, and R. Inger. 2018. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 6:e4794–32.

Sites web

GLMM for ecologists: Un site web très complet et avec une section Questions-Réponses pertinente !