Differences

This shows you the differences between two versions of the page.

--- r_atelier6 [2016/02/17 17:21]
zofia.taranu [GLM avec des données d'abondance]
+++ r_atelier6 [2016/09/05 16:44]
vincent_fugere
@@ Line 1: / Line 1: @@
-~~NOCACHE~~
+======= QCBS R Workshops =======
-======= Ateliers R du CSBQ =======
+[[http://qcbs.ca/|{{:logo_text.png?nolink&500|}}]]
-[[http://qcbs.ca/fr/|{{:logo_text.png?nolink&500|}}]]
 Cette série de [[r|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.
+====== Atelier 6: Modèles linéaires à effets mixtes ======
-====== Atelier 6: Modèles linéaires généralisés (mixtes) ======
+Développé par : Catherine Baltazar, Dalal Hanna, Jacob Ziegler
-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.
-**Résumé:** Les modèles linéaires généralisés sont des outils importants afin de surmonter un problème récurrent des modèles linéaires, c.-à-d. les variables de réponse n’ayant pas une distribution normale des résidus. Dans cet atelier, vous apprendrez les distributions principales en fonction de la nature de la variable de réponse, le concept de fonction de lien, et comment vérifier les suppositions de base de ces modèles. 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 présentation Prezi : [[https://prezi.com/vscqc2kaxnhu/|Prezi]]
-Lien vers la présentation Prezi associée : [[https://prezi.com/ycczxaot6tre/|Prezi]]
+Télécharger le script de R et les données pour cette leçon:
+  -[[http://qcbs.ca/wiki/_media/qcbs_w6_data.csv| Données des poissons]]
+  -[[http://qcbs.ca/wiki/_media/lmm_f.r| Code R]]
+===== Objectifs =====
+  - Qu’est-ce qu’un modèle linéaire à effets mixtes (MLM) et pourquoi est-ce important?
+  - Comment appliquer des MLMs dans R?
+     - Construction et exploration a priori des modèles et données
+     - Codage et sélection des modèles potentiels
+     - Validation des modèles
+     - Interprétation des résultats et visualisation du modèle
-Téléchargez le script R et les données pour cet atelier :
+===== Aperçu =====
-  - [[http://qcbs.ca/wiki/_media/script_atelier6.r| Script]]
-  - [[http://qcbs.ca/wiki/_media/mites.csv| Mites]]
-  - [[http://qcbs.ca/wiki/_media/faramea.csv| Faramea]]
-  - [[http://qcbs.ca/wiki/_media/banta_totalfruits.csv| Arabidopsis]]
-  - [[http://qcbs.ca/wiki/_media/inverts.csv| Invertébrés]]
-===== Limites des modèles linéaires et des modèles linéaires mixtes =====
-Pour illustrer l'utilité des modèles linéaires généralisés (GLMs en anglais), il est important de d'abord comprendre les limites des modèles linéaires (non-généralisés!), soit ceux présentés lors des ateliers 3 et 5. À cet effet, chargeons des données et appliquons quelques modèles linéaires:
+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.
-<file rsplus | charger les données>
+===== 1. Qu’est-ce qu’un MLM et pourquoi est-ce important? =====
-setwd("~/Desktop")
-mites <- read.csv('mites.csv')
-</file>
-Le jeu de données que vous venez de charger contient une partie du jeu de données des "mites Oribatidae". Ce dernier a été utilisé par plusieurs ouvrages sur R (p.ex. Borcard, Gillet & Legendre, //Numerical Ecology with R//) et est disponible dans la bibliothèque "vegan". Le jeu de données contient 70 échantillons de mousse et de mites provenant de la Station de Biologie des Laurentides de l'Université de Montréal, à Saint-Hippolyte, QC. Chaque échantillon contient des valeurs pour 5 variables environnementales de même que l'abondance de 35 taxa de mites. Le jeu de données que nous utiliserons pour cet atelier contient l'abondance d'un seul taxon de mite, "//Galumna sp.//", de même que les valeurs pour les 5 variables environnementales. Notre objectif est de modéliser l'abondance, la présence, et la proportion de Galumna en fonction des 5 variables environnentales. Nous avons donc créé une variable binaire présence/absence (1=présent, 0=absent) et une variable proportion (abondance de Galumna/abondance totale de mites).
+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.
-<file rsplus | explorer les données>
+====1.1 Introduction au jeu de données ====
-head(mites)
-str(mites)
-# 70 communautés de mites échantillonées à partir de carottes de mousse prises à la Station de Biologie des Laurentides, QC.
+<code rsplus | Chargez vos données >
-# Pour chaque carotte/échantillon, les valeurs suivantes sont fournies
+# Supprimer commandes antérieures en R
-# $Galumna: abondance de mites du genre Galumna
+rm(list=ls())
-# $pa: présence (1) ou absence (0) de Galumna, peu importe l'abondance
-# $totalabund: abondance totale de mites, toutes espèces confondues
-# $prop: proportion de Galumna dans la communauté de mites, i.e. Galumna/totalabund
-# $SubsDens: densité du substrat
-# $WatrCont: contenu d'eau de la carotte
-# $Substrate: type de substrat
-# $Shrub: abondance de buissons. Traité comme un facteur (une variable qualitative).
-# $Topo: microtopographie du sol. "couverture" (blanket) ou "hammac" (hummock).
-</file>
-Regardons si nous pouvons voir une relation entre Galumna et les variables environnementales à l'aide d'une simple matrice de diagrammes:
+# 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()
-<file rsplus | matrice de diagrammes>
+# Réglez le répertoire de travail vers le dossier qui contient les fichiers en copiant toute
-plot(mites)
+# la sortie de la commande file.choose (), à l'exception du nom de fichier R, et collez le tout dans setwd().
-</file>
-{{:6_spm.jpg?800|}}
+# 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()
-Il semble y avoir une relation négative entre l'abondance de Galumna et le contenu d'eau (WatrCont). La présence (pa) et la proportion (prop) de Galumna semblent aussi montrer une corrélation négative avec le contenu d'eau. Nous pouvons regarder ça de plus près en représentant dans des diagrammes les relations entre ces 3 variables réponses et avec le contenu d'eau:
+# 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)
-<file rsplus | 3 variables réponses vs. contenu d'eau>
+data <- read.csv('QCBS_W6_Data.csv')
-par(mfrow=c(1,3)) #division de la fenêtre de diagramme en une ligne et 3 colonnes pour avoir 3 diagrammes sur la même figure.
+</code>
-plot(Galumna ~ WatrCont, data = mites, xlab = 'Water content', ylab='Abundance')
-boxplot(WatrCont ~ pa, data = mites, xlab='Presence/Absence', ylab = 'Water content')
-plot(prop ~ WatrCont, data = mites, xlab = 'Water content', ylab='Proportion')
-par(mfrow=c(1,1)) #retour aux paramètres par défaut (un seul diagramme par figure)
-</file>
-{{:6_3plot.jpeg?600|}}
+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.
+{{:fig_1_qcbs_wiki.png?800|}}
-En effet, Galumna montre une relation négative avec le contenu d'eau, ce qui suggère que Galumna préfère les sites moins humides. Nous pouvons utiliser des modèles linéaires (fonction "lm") pour voir si ces relations sont statistiquement significatives.
+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.
-<file rsplus | modèles linéaires>
+----
-lm.abund <- lm(Galumna ~ WatrCont, data = mites)
+**DÉFI 1**
-summary(lm.abund)
-lm.pa <- lm(pa ~ WatrCont, data = mites)
-summary(lm.pa)
-lm.prop <- lm(prop ~ WatrCont, data = mites)
-summary(lm.prop)
-</file>
-Oui, il y a une forte relation significative pour les 3 variables réponses! Ça y est, on soummet à Nature. Attends une minute... Validons d'abord ces modèles pour voir s'ils respectent les suppositions des modèles linéaires, en commençant par le modèle d'abondance.
+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 :
+  - 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?
+  - S'attend-on à ce que ces relations soient pareilles entre les lacs ? Comment pourraient-elles différer ?
-<file rsplus | modèle d'abondance>
+++++ Réponse au défi 1 |
-plot(Galumna ~ WatrCont, data = mites)
+<code rsplus | Visualisez vos données >
-abline(lm.abund)
+# Utilisé pour simplifier les figures
-</file>
+fig <- theme_bw() + theme(panel.grid.minor=element_blank(), panel.grid.major=element_blank(),
+  panel.background=element_blank()) +
+  theme(strip.background=element_blank(), strip.text.y = element_text()) +
+  theme(legend.background=element_blank()) +
+  theme(legend.key=element_blank()) + theme(panel.border = element_rect(colour="black", fill=NA))
-{{:6_lm1.jpeg|}}
+# Faites les trois graphiques suivants pour explorer les données
+plot <- ggplot(aes(Fish_Length,Trophic_Pos),data=data)
-Le modèle ne représente pas très bien les données. Il prédit des valeurs négatives d'abondance lorsque le contenu d'eau excède 600, ce qui est insensé. En plus, le modèle ne prédit pas bien les valeurs élevées d'abondance lorsque le contenu d'eau est faible.
+# Graphique 1 - Toutes les données
+plot + geom_point() + xlab("Length (mm)") + ylab("Trophic Position") + labs(title="All Data") + fig
-<file rsplus | diagrammes diagnostiques>
+# Graphique 2 - Par espèce
-plot(lm.abund)
+plot + geom_point() + facet_wrap(~ Fish_Species) + xlab("Length (mm)") + ylab("Trophic Position") +
-</file>
+   labs(title="By Species") + fig
-{{:6_diagplots.jpeg|}}
+# Graphique 3 – Par lac
+plot + geom_point() + facet_wrap(~ Lake) + xlab("Length (mm)") + ylab("Trophic Position") +
+   labs(title="By Lake") + fig
+</code>
-Les diagrammes diagnostiques montrent que le modèle ne respecte pas les suppositions d'homogénéité de la variance (le diagramme à gauche montre que les rédidus sont plus grands pour les valeurs prédites élevées) ni de normalité (le diagramme à droite indique que les résidus ne sont pas distribués tel que prédit par une distribution normale, i.e. plusieurs points sont loin de la ligne pointillée). En conséquence, il nous faut rejeter ce modèle; nous ne pouvons pas l'utiliser pour conclure que l'abondance de Galumna varie en fonction du contenu d'eau. Les diagrammes diagnostiques indiquent que les modèles pour présence-absence et proportion sont eux aussi innapropriés:
+** SORTIE **\\
+** Graphique 1**
+{{:fig_2_w5.png?300|}}
+** Graphique 2**
+{{:fig_3_w5.png?300|}}
-<file rsplus | autres MLs>
-#Proportion
-plot(prop ~ WatrCont, data = mites)
-abline(lm.prop)
-plot(lm.prop)
-#Présence/Absence
-plot(pa ~ WatrCont, data = mites)
-abline(lm.pa)
-plot(lm.pa)
-</file>
-C'est très commun que les jeux de données biologiques ne respectent pas les suppositions d'homogénéité de la variance ou de normalité. Ces deux suppositions constituent les problèmes principaux des modèles linéaires, et sont la raison principale pourquoi nous avons besoin des modèles linéaires généralisés. Revoyons l'équation de base d'un modèle linéaire pour mieux comprendre d'où viennent ces suppositions. Cette équation est:
+** Graphique 3**
+{{:plot3.png?300|}}
-y<sub>i</sub> = β<sub>0</sub> + β<sub>1</sub>x<sub>i</sub> + ε<sub>i</sub>, où:
+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? Toutes les espèces semblent augmenter en position trophique avec la longueur, mais la pente peut être différente selon les espèces.
-  * y<sub>i</sub> est la valeur prédite pour la variable réponse,
+S'attend-on à ce que ces relations soient pareilles entre les lacs? Comment pourraient-elles différer? Certains paramètres spécifiques à un lac en particulier peuvent modifier la relation, tels que la productivité primaire du système.
-  * β<sub>0</sub> est l'ordonnée à l'origine de la droite de régression entre y et x,
-  * β<sub>1</sub> est la pente de la droite de régression entre y et x,
-  * x<sub>i</sub> est la valeur de la variable indépendante,
-  * ε<sub>i</sub> sont les résidus du modèle, qui proviennent d'une distribution normale avec une moyenne variable mais une variance constante.
-Ce dernier point à propos de ε<sub>i</sub> est important. C'est de là que proviennent les suppositions de normalité et d'homogénéité de la variance (homoscédasticité). Ceci veut dire que les residus du modèle (les distances entre chaque observation et la droite de régression) peuvent être prédits en pigeant de façon aléatoire des valeurs dans une distribution normale. Souvenez-vous que toutes les distributions/lois normales ont deux paramètres, soit μ (la moyenne de la distribution) et σ<sup>2</sup> (la variance de la distribution):
+++++
-{{:6_normal_params.jpeg|}}
-Dans un modèle linéaire, μ change selon la valeur de x (la variable indépendante), mais σ<sup>2</sup> garde la même valeur pour toutes les valeurs de y<sub>i</sub>. En effet, une autre équation pour représenter un modèle linéaire est la suivante:
-y<sub>i</sub> ~ //N//(μ = β<sub>0</sub> + β<sub>1</sub>x<sub>i</sub>, σ<sup>2</sup>),
+====1.2 Comment pourrions-nous analyser ces données?====
-ce qui signifie littéralement que chaque observation (y<sub>i</sub>) provient d'une distribution normale avec les paramètres μ (qui dépend de la valeur de x<sub>i</sub>) et σ<sup>2</sup>. Éssayons de prédire l'abondance de Galumna en fonction du contenu d'eau en utilisant cette équation et les paramètres que nous avons estimé plus tôt avec la fonction lm(). D'abord, il nous faut des valeurs pour β<sub>0</sub> et β<sub>1</sub> (les coéfficients de régression), que nous pouvons obtenir ainsi:
+=== 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 :
-<file rsplus | coeffs ML>
+{{:fig_5_w5.png?400|}}
-coef(lm.abund)
-#(Intercept)     WatrCont
-#3.439348672 -0.006044788
-</file>
-Ce modèle prédit que pour un contenu d'eau de disons 300, nous devrions obtenir une abondance de Galumna de 1.63:
+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.
-.44 + (-0.006 x 300) = 1.63.
+=== 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 :
-.63 est l'abondance prédite s'il n'y avait pas de variance résiduelle (i.e. si notre modèle avait un r2 de 1). Pour modéliser les valeurs prédites, nous devons donc ajouter ε<sub>i</sub>, i.e. la distance entre les valeurs observées et la droite de régression. C'est ici que nous utilisons la distribution normale. Pour x = 300, notre modèle prédit que ε<sub>i</sub> devrait suivre une distribution normale avec une moyenne = 1.63. Lorsque le contenu d'eau = 400, ε<sub>i</sub> devrait suivre une distribution normale avec une moyenne =  3.44 + (-0.006 x 400) = 1.02. Chaque valeur de y<sub>i</sub> est modélisée en utilisant une distribution normale différente, chacune avec une moyenne distincte qui dépend de x<sub>i</sub>. Par contre, la variance de toutes ces distributions (σ<sup>2</sup>) reste toujours la même. La fonction lm() trouve la valeur de σ<sup>2</sup> optimale pour minimiser la somme des carrés résiduelle et utilise ensuite cette valeur pour toutes les distributions normales servant à modéliser y. Nous pouvons trouver la valeur de σ<sup>2</sup> de notre modèle dans le résumé du modèle:
+{{:fig_6_w5.png?400|}}
-<file rsplus | sigma>
+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.
-summary(lm.abund)$sigma
-</file>
-Nous trouvons que sigma est plus ou moins 1.51. Nous avons maintenant tous les coéfficients nécessaires pour modéliser manuellement l'abondance de Galuman comme le ferait la fonction lm(). À un contenu d'eau de 300, les résidus devraient suivre une distribution normale avec les paramètres μ = 1.63 et σ<sup>2</sup> = 1.51. À un contenu d'eau de 400, les résidus devraient suivre une distribution normale avec les paramètres μ = 1.02 et σ<sup>2</sup> = 1.51, etc. Graphiquement, ça ressemble à ça:
+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).
-{{:6_lm_assump.jpeg|}}
+===Construire un MLM===
+Les MLMs sont un compromis entre séparer et regrouper. Ils:
+  - 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.
+  - Utilisent toutes les données disponibles (regrouper) tout en contrôlant les différences entre les lacs et les espèces (pseudo-réplication).
-Les quatres lois normales sur ce graphique donnent la probabilité d'observer toutes valeurs possibles de Galumna à 4 contenus d'eau différents. La moyenne de la distribution varie selon le contenu d'eau (donc μ diminue avec le contenu d'eau), mais σ<sup>2</sup> est toujours 1.51. La **variance est donc homogène** pour toutes les valeurs de x. Ce modèle est inadéquat pour au moins deux raisons:
+====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.
-. Les valeurs prédites sont en moyenne plus loin de la droite de régression lorsque le contenu d'eau est faible. Ceci signifie que la variance résiduelle est plus grande lorsque x est faible, ce qui revient à dire que ε<sub>i</sub> varie en fonction de x et que la variance n'est donc pas homogène. Il est innaproprié ici d'utiliser une valeur constante de σ<sup>2</sup>: les distributions normales utilisées pour prédire y devraient être plus larges (avoir une plus grande variance) lorsque x est faible que lorsqu'il est élevé. Malheureusement, les modèles linéaires ne permettent pas ceci.
+=== 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.
-. C'est innaproprié d'utiliser une distribution normale pour prédire y en fonction de x. Notre variable réponse est l'abondance qui, par définition, doit être un nombre entier. Pourtant, lorsque le contenu d'eau = 300, notre modèle prédit que l'abondance la plus probable est 1.63! Nous savons pourtant que la probabilité d'observer une telle valeur (ou toute autre fraction) est 0. Nos valeurs prédites devraient être modélisées en utilisant une distribution qui ne contient que des nombres entiers, plutôt qu'avec une distribution continue telle que la distribution normale. Ceci est un problème courant car les variables biologiques peuvent être distribuées de bien d'autres façon que ce qui est prédit par la distribution normale.
+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.
-Les modèles linéaires généralisés peuvent résoudrent ces deux problèmes. Poursuivez votre lecture!
+=== Effet aléatoire===
-===== Les données biologiques et leurs distributions =====
+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.
-Les statisticiens ont développé une foule de lois de probabilité (ou distributions) correspondant à divers types de données. Une **loi de probabilité** donne la **probabilité** d'observer chaque **issu** possible d'une expérience ou campagne d'échantillonage (par. ex. abondance = 8 Galumna est un issu d'un échantillonage). Les lois “discrètes” n'incluent que des nombres entiers dans leur ensemble d'issus, alors que les lois “continues” incluent aussi des fractions (par. ex. la loi normale). Toutes les lois ont des paramètres qui déterminent la forme de la loi/distribution (par. ex. μ et σ2 pour la loi normale). Pour un excellent survol des lois de probabilités utiles en écologie, nous vous recommandons le chapitre 4 du livre de Ben Bolker //Ecological Models and Data in R//. Ici nous ne présenterons que brièvement quelques distributons utiles pour les GLMs.
+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.
-Nous avons déjà vu que notre variable “abondance de Galumna ” ne peut prendre comme valeur que des nombres entiers, et serait donc mieux modélisée par une loi discrète qu'une loi continue. Une loi utile pour modéliser les données d'abondance est la loi de "Poisson", nommé en l'honneur du statisticien Siméon Denis Poisson. La loi de Poisson est une loi discrète avec un seul paramètre:  λ (lambda), qui détermine et la moyenne et la variance de la distribution (la moyenne et la variance d'une loi de Poisson sont donc égales). Voici 3 exemples de lois de Poisson avec des valeurs différentes de λ, ce qui correspond ici au nombre moyen de Galumna retrouvé dans un ensemble fictif d'échantillons:
-{{:6_poisson.jpeg|}}
+====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).
-Remarquez que lorsque λ est faible (c.-à-d. lorsque la moyenne s'approche de zéro), la distribution est décalée vers la gauche, alors que lorsque λ est élevée, la distribution est symmétrique. La variance augmente aussi avec la moyenne (puisque les deux ont la même valeur), les valeurs prédites sont toujours des nombres entiers, et l'ensemble d'issus d'une loi de Poisson est strictement positif. Toutes ces propriétés sont utiles pour modéliser les données de dénombrement, p.ex. l'abondance d'un taxon, le nombre de graines dans une parcelle, etc. Notre variable mites$Galumna semble suivre une loi de Poisson avec une basse valeur de λ (en effet, si nous calculons l'abondance moyenne de Galumna pour l'ensemble des échantillons avec la fonction mean(), nous observons que cette valeur est proche de 0):
+**Intercepts**
-<file rsplus | histogramme pour abondance de Galumna>
+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).
-hist(mites$Galumna)
-mean(mites$Galumna)
-</file>
-{{:6_galumnahist.jpeg|}}
+{{:fig_7_w5.png?600|}}
-La variable mites$pa (présence-absence) prend une autre forme. Cette variable n'inclue que des 0s et des 1s, de telle sorte que la loi de Poisson ne serait pas plus appropriée pour cette variable que la loi normale.
+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).
-<file rsplus | histogramme pour présence-absence de Galumna>
+{{:fig_8_w5.png?600|}}
-hist(mites$pa)
-</file>
-{{:6_pa.jpeg|}}
+**Pentes**
-Nous avons besoin d'une distribution qui n'inclut dans son ensemble que deux issus possibles: 0 ou 1. La loi de “Bernoulli” est une distribution de la sorte. C'est souvent la première loi de probabilité qu'on nous présente dans les cours de statistiques, pour prédire la probabilité d'obtenir "pile" ou "face" en tirant à... pile ou face. Cette loi n'a qu'un paramètre: //p//, la probabilité de succès. Nous pouvons utiliser la loi de Bernouilli pour calculer la probabilité d'obtenir l'issu “Galumna présent” (1) vs. “Galumna absent” (0). Voici des exemples de distributions de Bernoulli avec différentes probabilités de présence (//p//):
+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.
-{{:6_bernouill.jpeg|}}
+{{:fig_9_w5.png?600|}}
-Nous pouvons calculer le nombre de sites où Galumna est présent par rapport au nombre total de sites pour obtenir un estimé de ce que //p// pourrait être dans cet exemple:
+===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.
-<file rsplus | p pour Galumna présence-absence>
+{{:fig_12_w5.png?600|}}
-sum(mites$pa) / nrow(mites)
-</file>
-//p// pour la variable mites$pa est plus ou moins 0.36, de tel sorte qu'environ deux fois plus de sites montrent l'issu "Galumna absent" (0) que l'issu "Galumna présent" (1).
+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.
-Lorsqu'il y a plusieurs épreuves (chacune avec un succès/échec), la loi de Bernoulli se transforme en loi binomiale, qui inclue le paramètre additionel //n//, le nombre d'épeuves. La loi binomiale prédit la probabilité d'observer une certaine proportion de succès, //p//, sur le nombre total d'épreuves, //n//. Un "succès" pourrait être, par exemple, la présence d'un taxon à un site (comme pour Galumna), le nombre d'individus qui survit pendant une année, etc. Voici des exemples de lois binomiales avec //n// = 50 et trois valeurs différentes de //p//:
+**CIC élevé**
-{{:6_binomial.jpeg|}}
+{{:fig_10_w5.png?400|}}
-Remarquez que la loi binomiale est assymétrique et décalée à gauche lorsque //p// est faible, mais elle est décalée à droite lorsque //p// est élevé. C'est la différence principale avec la loi de Poisson: l'étendu de la loi binomial a une limite supérieure, correspondant à //n//. Conséquemment, la loi binomiale est utilisée pour modéliser des données lorsque le nombre de succès est donné par un nombre entier, et lorsque le nombre d'épreuves est connu. Nous pouvons utiliser la distribution binomiale pour modéliser nos données de proportion, où chaque mite échantillonée pourrait être considérée comme une épreuve. Si la mite est du genre Galumna, l'épreuve est un succès (1), sinon, c'est un échec (0). Dans ce cas, le nombre d'épreuves //n// varie pour nos 70 échantillons selon l'abondance totale de mites dans l'échantillon, alors que //p//, la probabilité de succès, est donné par la proportion de Galumna dans chaque échantillon.
+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.
-Pourquoi tout cette discussion à propos des lois de distribution? Parce que n'importe quelle loi peut être utilisée pour remplacer la loi normale lorsqu'on calcule les valeurs estimées dans un modèle linéaire. Par exemple, nous pouvons utiliser la loi de Poisson pour modéliser nos valeurs d'abondance en utilisant l'équation suivante:
+**CIC faible**
-y<sub>i</sub> ~ Poisson(λ = β<sub>0</sub> + β<sub>1</sub>x<sub>i</sub>)
+{{:fig_11_w5.png?400|}}
-Remarquez que λ varie selon x (contenu d'eau), ce qui signifie que la variance dans les résidus variera aussi selon x (car pour la loi de Poisson variance = aussi λ). Ceci veut dire que nous venons de nous défaire de la supposition d'homogénéité de la variance! Aussi, les valeurs estimées seront désormais des nombres entiers plutôt que des nombres décimaux car ils seront tirés de lois de Poisson avec différentes valeurs de λ. Ce modèle ne prédiera jamais de valeurs négatives car l'ensemble d'une loi de Poisson est toujours strictement positif. En changeant la distribution des résidus (ε<sub>i</sub>) de normale à Poisson, nous avons corrigé presque tous les problèmes de notre modèle linéaire pour l'abondance de Galumna.
+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.
-Ce modèle est __presqu__'un GLM de Poisson, qui ressemble à ça:
+----
+**DÉFI 2**
-{{:6_poissonglm.jpeg|}}
+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?
+  - Les positions trophiques des poissons ne varient pas entre les lacs?
+  - Les positions trophiques des poissons sont similaires dans les lacs mais différentes entre les lacs?
-Remarquez que les probabilités d'observer différentes valeurs estimées (en orange) sont maintenant des nombres entiers, et qu'autant la variance que la moyenne de la distribution declinent lorsque λ diminue avec le contenu d'eau. Pourquoi la droite de régression est-elle courbée? Pourquoi est-ce que ce modèle se nomme un "modèle linéaire généralisé"? Continuez votre lecture!
+++++ Réponse défi 2 |
+  - CIC faible - intervalles de confiance plus petits
+  - CIC élevé - intervalles de confiance plus larges
+++++
+----
+=====2. Comment implémenter un MLM dans R?=====
+Le protocole des modèles mixtes dans R:
+.1: Construction du modèle a priori et exploration des données
+.2: Coder les modèles potentiels et sélectionner le meilleur modèle
+.3: Valider le modèle
+.4: Interpréter les résultats et visualiser le modèle
-===== GLM avec une variable réponse binaire =====
+====2.1 Construction du modèle a priori et exploration des données====
-Les variables binaires sont fréquentes en écologie : on observe un phénomène X ou son "absence". Par exemple, on note souvent la présence ou l'absence d'espèces lors d'études de terrain. Le but est généralement de déterminer si la présence d'une espèce est influencée par différentes variables environnementales. D'autres exemples courants sont la présence/absence d'une maladie au sein d'une population sauvage, l'observation/non-observation d'un comportement spécifique et la survie/mort d'un individu. Un modèle de régression qui utilise une variable binaire comme variable réponse est l'un de plusieurs modèles linéaires généralisés (GLM) et est appelé régression logistique ou modèle logit.
+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:
-Dans R, la présence (ou succès, survie…) est habituellement codée par un ''1'' et une absence (ou échec, mort…) par un ''0''. On effectue une régression logistique (ou n'importe quel GLM) à l'aide de la fonction ''glm()''. Cette fonction diffère un peu de la fonction de base ''lm()'', car elle permet de spécifier une distribution statistique autre que la distribution normale. Nous avons déjà vu que les variables binaires ne sont pas distribuées normalement (//i.e.// on observe un pic à 0 et un pic à 1 et rien entre les deux). La dernière section a montré que la distribution de Bernoulli est appropriée pour modéliser les variables réponses binaires. La moyenne de cette distribution est représentée par la probabilité //p// d'observer un résultat et la variance est calculée par //p//*(1 - //p//). Le terme (1 - //p//) représente la probabilité de ne **pas** observer un résultat. Dans R, on spécifie la distribution statistique du modèle avec l'argument ''family''. Pour la régression logistique, on l'indique de la façon suivante : ''family = 'binomial'''. Rappelez-vous que la distribution de Bernoulli est un cas spécifique de la distribution binomiale lorsque le nombre de répétitions est égal à 1: R "comprend" qu'il faut utiliser une distribution de Bernoulli.
+<m> {PT_ijk} ∼ {Longueur_i} + {Lac_j} + {Espèce_k} + {ε} </m>
-Lorsqu'on prédit la probabilité d'observer un phénomène Y qui est une variable binaire, la valeur prédite doit se trouver entre 0 et 1 : c'est l'étendue possible d'une probabilité ! Si on utilise un modèle linéaire de base pour modéliser une variable binaire en fonction de variables explicatives, il est fort possible qu'on obtienne des valeurs prédites qui se retrouvent en dehors de l'intervalle [0,1], ce qui ne fait aucun sens. L'exemple suivant va vous permettre de mieux comprendre pourquoi un modèle linéaire traditionnel est inapproprié. La prochaine sous-section va vous montrer comment éviter ce problème avec une fonction de lien. Brièvement, une fonction de lien est utilisée pour linéariser la relation entre les valeurs prédites d'un modèle et le prédicteur linéaire (voir sous-section suivante).
+où,\\
-<file rsplus | Utilisation inappropriée d'un modèle linéaire>
+<m>{PT_ijk}</m> est la position trophique du poisson (i) du lac (j) et de l’espèce (k)\\
-model.lm <- lm(pa ~ WatrCont + Topo, data = mites)
+et\\
-fitted(model.lm)
+ε sont les résidus du modèle (c. à d. la variation inexpliquée).\\
-# La fonction "fitted()" extraie les valeurs prédites de la variable réponse du modèle linéaire.
-# Certaines valeurs sont en-dessous de 0, ce qui ne fait pas de sens pour une régression logistique.
-# Essayong le même modèle, mais avec une distribution binomiale cette fois-ci.
-# Remarquez l'argument "family" pour spécifier la distribution.
-model.glm <- glm(pa ~ WatrCont + Topo, data = mites, family = binomial)
-fitted(model.glm)
-# Toutes les valeurs sont comprises entre 0 et 1.
-</file>
-==== Le concept de la fonction de lien ====
+**Exploration des données **
-Afin d'éviter les biais reliés aux modèles linéaires de base, nous avons besoin de spécifier deux choses : une distribution statistique pour les résidus du modèle et une fonction de lien pour les valeurs prédites par ce même modèle. Nous avons déjà vu la distribution de Bernoulli dans la section précédente, alors nous passerons directement à l'explication du rôle de la fonction de lien.
+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.
-Pour un modèle de régression linéaire d'une variable réponse continue distribuée normalement, l'équation suivante nous permet d'obtenir les valeurs prédites de la variable réponse :
+<code rsplus | 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
-//μ// = //Xβ//
+# 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
-où //μ// est la valeur prédite de la réponse variable, //X// est la matrice du modèle (//i.e.// ça représente les variables explicatives) et //β// correspond aux paramètres estimés à partir des données (//i.e.// l'intercept et la pente). Le terme //Xβ// est appelé le prédicteur linéaire. En termes mathématiques, c'est le produit matriciel de la matrice du modèle //X// et du vecteur des paramètres estimés //β//. Regardons cela de plus près dans R :
+# ii)Housekeeping et exploration des données
+# Assurez-vous que la structure de vos données soit correcte
+str(data)
+</code>
-<file rsplus | Illustration du concept de prédicteur linéaire>
+++++ Sortie |
-# Chargez le jeu de données CO2. C'est ce jeu de données que nous avons utilisé dans l'atelier 3 !
+{{:fig_13_w5.png?600|}}
-data(CO2)
+++++
-head(CO2)
-# Construisez un modèle linéaire de l'absorption de CO2 en fonction de la concentration ambiante de CO2.
-model.CO2 <- lm(uptake ~ conc, data = CO2)
-# On extraie la matrice du modèle avec la fonction model.matrix().
-X <- model.matrix(model.CO2)
-# Les paramètres estimés sont extraits ainsi :
-B <- model.CO2$coefficients
-# On multiple X et B pour obtenir le prédicteur linéaire.
-# Le symbole "%*%" indique qu'on veut effectuer le produit matriciel.
-XB <- X %*% B
-# On compare les valeurs de XB aux valeurs obtenues avec la fonction fitted().
-# Toutes les déclarations devraient être vraies (i.e. TRUE).
-# On utilise la fonction round() pour que tous les éléments aient cinq décimales.
-round(fitted(model.CO2), digits = 5) == round(XB, digits = 5)
-</file>
-Lorsqu'on crée un modèle linéaire simple avec une variable réponse continue distribuée normalement, le prédicteur linéaire est égal aux valeurs attendues de la variable réponse. Ceci n'est pas exact si la variable réponse n'est pas distribuée normalement. Si c'est le cas, il faut appliquer une transformation sur les valeurs prédites, //i.e.// une fonction de lien. La fonction de lien peut être vue comme une transformation des valeurs prédites pour obtenir une **relation linéaire** entre celles-ci et le prédicteur linéaire :
+Regardez la distribution des échantillons pour chaque facteur :
+<code rsplus | 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)
+</code>
-//g//(//μ//) = //Xβ//
+++++ Sortie |
+**Lac**\\
+----
+{{:fig_15_w5.png?200|}}\\
+----
+**Espèce**\\
+----
+{{:fig_14_w5.png?100|}}\\
+----
+++++
-où //g//(//μ//) est la fonction de lien des valeurs prédites. Ceci permet d'enlever la contrainte de distribution normale des résidus. Lorsque la variable réponse est une variable binaire, la fonction de lien est la fonction logit et est représentée par l'équation suivante :
+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!)
-logit(//μ//) = log (//μ// / 1-//μ//) = //Xβ//
+Regardez la distribution des variables continues :
+<code rsplus | 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)
+</code>
-où //μ// représente les valeurs prédites (//i.e.// la probabilité que Y = 1, car on observe la présence d'une espèce, d'une maladie, d'un succès ou d'un autre événement). Le ratio //μ// / 1-//μ// représente la cote (odds en anglais) qu'un événement se produise. Cela transforme les valeurs prédite sur une échelle de 0 à l'infini. Par exemple, si on a une probabilité de 0,8 d'observer une espèce X, la cote est donc de 4 : il est 4 fois plus probable d'observer l'espèce que de ne pas l'observer - 0.8/(1-0.8) = 4. La transformation log (on appelle maintenant ce ratio le log odds) permet aux valeurs d'être distribuées de moins l'infini à l'infini. Donc, la fonction de lien logit prend les valeurs prédites du modèle et les transforme en une variable continue sans borne. Les valeurs prédites peuvent ainsi être reliées à un prédicteur linéaire. C'est pour cette raison qu'on appelle ce modèle un modèle **linéaire** généralisé même si la relation entre la variable réponse et la variable explicative ne ressemble pas nécessairement à une "ligne droite" !
+++++ Sortie |
+**Longueur**\\
-<file rsplus | Un exemple simple de régression logistique>
+----
-# On construit un modèle de régression de la présence/absence d'une espèce de mite (Galumna sp.)
+{{:fig_16_w5.png?400|}}\\
-# en fonction du contenu en eau du sol et de la topographie.
+----
-# On utilise la fonction glm() et on spécifie l'argument "family".
+**Position trophique**\\
-logit.reg <- glm(pa ~ WatrCont + Topo, data = mites, family = binomial(link = "logit"))
+----
-# La fonction de lien "logit" est la fonction par défaut pour une régression logistique,
+{{:fig_17_w5.png?400|}}\\
-# ce qui signifie qu'il n'est pas nécessaire de l'indiquer avec l'argument "family":
+----
-logit.reg <- glm(pa ~ WatrCont + Topo, data = mites, family = binomial)
-summary(logit.reg)
-</file>
-** DÉFI 5 **
-En utilisant le jeu de données ''bacteria'' (du paquet ''MASS''), modélisez la présence de //H. influenzae// en fonction du traitement (''trt'') et de la semaine de test (''week''). Commencez avec un modèle saturé et réduisez-le jusqu'à obtenir le modèle le plus parcimonieux possible.
-++++ Défi 5 : Solution |
-<file rsplus | >
-model.bact1 <- glm(y ~ trt * week, family = binomial('logit'), data = bacteria)
-model.bact2 <- glm(y ~ trt + week, family = binomial('logit'), data = bacteria)
-model.bact3 <- glm(y ~ week, family = binomial('logit'), data = bacteria)
-anova(model.bact1, model.bact2, model.bact3, test = 'LRT')
-# Analysis of Deviance Table
-# Model 1: y ~ trt * week
-# Model 2: y ~ trt + week
-# Model 3: y ~ week
-#   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-# 1       214     203.12
-# 2       216     203.81 -2  -0.6854  0.70984
-# 3       218     210.91 -2  -7.1026  0.02869 *
-# En se basant sur ces résultats, on doit choisir le modèle 2
-# comme le modèle représentant le mieux le jeu de données.
-</file>
 ++++
-==== Interpréter la sortie d'une régression logistique ====
+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.
-La sortie du modèle de régression logistique indique que les deux variables explicatives (''WatrCont'' et ''Topo'') sont significatives, mais comment interprète-on les coefficients estimés ? Rappelez-vous que nous avons effectué une transformation des valeurs prédites par le modèle (//i.e.// la probabilité que Y = 1), alors il faut utiliser une fonction inverse pour pouvoir interpréter correctement les résultats. On peut utiliser la fonction exponentielle ''e<sup>x</sup>'' pour obtenir la cote de probabilité pour chaque variable explicative.
+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.
-<file rsplus | Interpréter la sortie d'une régression logistique>
+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 :
-# Pour obtenir les pentes, il faut utiliser la fonction exponentielle exp().
-# Ceci mettra les coefficients sur l'échelle des cotes.
-# Mathématiquement, ceci correspond à :
-# exp(model coefficients) = exp(log(μ / (1 - μ)) = u / (1 - μ)
-# C'est notre rapport de cotes !
-exp(logit.reg$coefficients[2:3])
-#  WatrCont    TopoHummock
-#  0.9843118   8.0910340
-# Pour obtenir l'intervalle de confiance sur l'échelle des cotes :
-exp(confint(logit.reg)[2:3,])
-#               2.5 %      97.5 %
-#  WatrCont     0.9741887  0.9919435
-#  TopoHummock  2.0460547  38.6419693
-</file>
-Prenez note que la cote pour une variable explicative est calculée lorsque les autres variables sont gardées constantes. La topographie a une cote de 8.09. Ceci signifie que la probabilité d'observer //Galumna// sp. est 8.09 fois plus vraisemblable lorsque la topographie est de type ''hummock'' plutôt que ''blanket''.
+<code rsplus | Exploration de données >
+# Évaluer la colinéarité entre variables
+plot(data)
+cor(data$Fish_Length, data$Var2)
+</code>
-Lorsque la cote est inférieure à 1, l'interprétation est un peu plus compliquée. Si c'est le case, il faut prendre la valeur inverse de la cote (//i.e.// 1 divisé par la cote) pour faciliter l'interprétation. L'interprétation revient à dire comment l'observation d'un phénomène est **MOINS** probable. Pour le contenu en eau du sol, la cote est de 0.984. L'inverse est 1 / 0.984 = 1.0159. Ceci signifie que l'augmentation d'une unité en contenu en eau diminue la vraisemblance d'observer la présence de //Galumna// sp. de 1.0159. On peut aussi l'exprimer en pourcentage en soustrayant 1 à cette valeur : (1.0159 - 1) * 100 = 1.59 %. Il est 1.59 % moins vraisemblable d'observer //Galumna// sp. avec une augmentation d'une unité de contenu en eau. Pour se convaincre qu'on fait la bonne interprétation, on peut représenter graphiquement les résultats de la présence de //Galumna// sp. en fonction du contenu en eau du sol. On voit qu'en moyenne la présence de //Galumna// sp. est plus élevée lorsque le contenu en eau est faible
+----
+**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 ?
-{{:galumna_pa.png?400|}}
+++++ Réponse défi 3 |
+Un exemple est la masse du poisson – c’est une variable fortement corrélée avec la longueur du poisson. Par conséquent, nous ne voulons pas inclure ces deux variables dans le même modèle.
+++++
+----
-Lorsqu'un paramètre estimé est entre 0 et 1 sur l'échelle des cotes, la relation entre la variable réponse et la variable explicative est négative. Si la valeur est supérieure à 1, cela indique une relation positive entre les variables. Si l'intervalle de confiance inclut la valeur 1, la relation entre les variables n'est pas significative. Rappelez-vous qu'une valeur de 1 sur l'échelle des cotes signifie que la probabilité d'observer un phénomène Y est la même que celle de ne pas observer ce phénomène (//i.e.// quand p = 0.5, 0.5/(1-0.5) = 1).
+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 :\\
+<m> {z} = ({x} - {mean(x)}) / {sd(x)} </m>\\
-Pour obtenir une probabilité au lieu d'une cote pour chaque variable explicative, il faut utiliser la fonction logit inverse :
+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.
+<code rsplus | 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)
+</code>
-logit<sup>-1</sup> = 1/(1+1/exp(x))
+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 :
+  - Créant un modèle linéaire sans les facteurs qui pourraient avoir un effet aléatoire
+  - Calculant les résidus de ce modèle linéaire
+  - Produisant un graphique de la valeur des résidus en fonction des niveaux des facteurs potentiellement aléatoires
-où x est le paramètre à transformer de l'échelle log odds à l'échelle de probabilité. Pour le modèle ''logit.reg'', le paramètre estimé pour la topographie est de 2.091 sur l'échelle log odds. Donc, la probabilité est donnée par :
+<code rsplus | 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)
+</code>
-/(1+1/exp(2.091)) = 0.89 ce qui équivaut à 1/(1+1/8.09). Rappelez-vous que la valeur 8.09 est sur l'échelle des cotes. On a une probabilité de 0.89 d'observer //Galumna// sp. lorsque la topographie est de type ''hummock''.
+++++ Sortie|
+**Effet de l’espèce**\\
+----
+{{:fig_18_w5.png?400|}}\\
+----
+**Effet du lac**\\
+----
+{{:fig_19_w5.png?400|}}\\
+----
+++++
-<file rsplus | Un peu de mathématique sur la fonction logit inverse>
+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.
-# On commence avec la valeur de cote pour la topographie du modèle logit.reg :
-µ/ (1 - µ) = 8.09
-# On réarrange pour isoler µ :
-µ = 8.09(1 - µ) = 8.09 - 8.09µ
-.09µ + µ = 8.09
-µ(8.09 + 1) = 8.09
-µ = 8.09 / (8.09 + 1)
-µ = 1 / (1 + (1 / 8.09)) = 0.89
-# On obtient le même résultat sans utiliser la fonction logit inverse !
-</file>
-==== Pouvoir prédictif et validation du modèle ====
-Une façon simple et intuitive d'estimer le pouvoir explicatif d'un GLM est de comparer la déviance du modèle à celle d'un modèle nul. La déviance peut être vue comme une généralisation du concept de la somme des carrés résiduelle lorsque le modèle est estimé par maximisation de la vraisemblance (//i.e.// la méthode par laquelle on estime les paramètres d'un GLM). Ceci nous permet de calculer un pseudo-R<sup>2</sup>, une statistique similaire au R<sup>2</sup> dans une régression des moindres carrés (//i.e.// la méthode utilisée pour estimer les paramètres d'une régression linéaire de base). Le modèle nul correspond à un modèle sans variable explicative. Dans R, on l'indique de la façon suivante : ''null.model %%<-%% glm(Response.variable ~ 1, family = binomial)''. La forme générique pour calculer un pseudo-R<sup>2</sup> est :
+====2.2 Coder les modèles potentiels et sélectionner le meilleur modèle ====
+Notre modèle a priori \\
-Pseudo-R<sup>2</sup> = (déviance du modèle nul – déviance résiduelle) / déviance résiduelle
+<m> {PT_ijk} ∼ {Longueur_i} + {Lac_j} + {Espèce_k} + {ε} </m>\\
-où "déviance du modèle nul" est la déviance du modèle nul et "déviance résiduelle" est la déviance résiduelle du modèle d'intérêt. La différence est divisée par la déviance du modèle nul afin de contraindre le pseudo-R<sup>2</sup> entre 0 et 1.
+Dans R, on le code ainsi :\\
+{{:fig_20_w5.png?500|}}
-<file rsplus | Le pseudo-R2 d'une régression logistique>
+La structure de lmer n’est pas intuitive. Les éléments de base de la fonction sont :\\
-# Les déviances résiduelle et nulle sont déjà enregistrées dans un objet de type glm.
+{{:fig_21_w5.png?500|}}\\
-objects(logit.reg)
-pseudoR2 <- (logit.reg$null.deviance – logit.reg$deviance) / logit.reg$null.deviance
-pseudoR2
-# [1]  0.4655937
-</file>
-Les variables explicatives du modèle expliquent 46.6% de la variabilité de la variable réponse.
+REML (Restricted Maximum Likelihood) est la méthode par défaut dans la fonction "lmer".
-Récemment, [[http://www.tandfonline.com/doi/abs/10.1198/tast.2009.08210#.VFpKZYcc4ow|Tjur (2009)]] a proposé une nouvelle statistique, le coefficient de discrimination (//D//), afin d'évaluer le pouvoir prédictif d'une régression logistique. Intuitivement, //D// évalue si la régression logistique peut classer adéquatement chaque résultat comme un succès ou un échec. Mathématiquement, c'est la différence entre la moyenne des valeurs prédites des succès (//i.e.// Y = 1) et des échecs (//i.e.// Y = 0) :
+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).
-<m>D = overline{π}_1 - overline{π}_0</m>
+Mais comment faire si on souhaite aussi que la pente puisse varier?\\
+{{:fig_22_w5.png?500|}}\\
-où <m>overline{π}_1</m> est la moyenne attendue des probabilités lorsqu'un événement est observé et <m>overline{π}_0</m> est la moyenne attendue des probabilités lorsqu'un événement n'est pas observé. Une valeur de //D// près de 1 indique que le modèle accorde une probabilité élevée d'observer un événement lorsque celui-ci a été observé et une probabilité faible d'observer un événement lorsque celui-ci n'a pas été observé. Une valeur de //D// près de 0 indique que le modèle n'est pas utile pour distinguer entre les observations et les "non-observations" d'un résultat. Le code suivant montre comment obtenir la valeur de //D// et comment représenter visuellement  les valeurs de <m>overline{π}_1</m> et <m>overline{π}_0</m>.
+----
+**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.
-<file rsplus | Le coefficient de discrimination et sa visualisation>
+<code rsplus | Coder les MLMs>
-install.packages("binomTools")
+lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)
-library("binomTools")
+</code>
-# La fonctions Rsq calcule indices d'ajustement
+++++ Réponse défi 4 |
-# dont le coefficient de discrimination.
+<code rsplus | Coder MLMs>
-# Pour plus d'information sur les autres indices, consultez Tjur (2009).
+lmer(Z_TP~Z_Length + (1+Z_Length|Lake) + (1+Z_Length|Species), data = data, REML=TRUE)
-# Les histogrammes montrent la distribution des valeurs attendues lorsque le résultat est observé
+</code>
-# et non observé.
+++++
-# Idéalement, le recouvrement entre les deux histogrammes doit être minimal.
+----
-fit <- Rsq(object = logit.reg)
-fit
-# R-square measures and the coefficient of discrimination, 'D':
-#
-#    R2mod     R2res     R2cor     D
-#    0.5205221 0.5024101 0.5025676 0.5114661
-#
-# Number of binomial observations:  70
-# Number of binary observation:  70
-# Average group size:  1
-plot(fit, which = "hist")
-</file>
-Pour évaluer l'ajustement (//i.e.// goodness-of-fit) d'une régression logistique, les graphiques de diagnostic ne sont pas très utiles (voyez l'atelier 3). On utilise plutôt un [[http://en.wikipedia.org/wiki/Hosmer-Lemeshow_test|test de Hosmer-Lemeshow]] pour évaluer si un modèle est approprié ou non. Ce test sépare les valeurs attendues (ordonnées de la plus faible à la plus grande) en groupes de même taille. Dix groupes est généralement le nombre de groupes recommandé. Pour chaque groupe, on compare les valeurs observées aux valeurs attendues. Le test est similaire à un test de khi-carré avec G - 2 degrés de liberté (G est le nombre de groupes). Dans R, ce test est disponible dans le paquet ''binomTools''.
+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».
-<file rsplus | Le test de Hosmer-Lemeshow>
+<code rsplus | Coder les MLMs>
-fit <- Rsq(object = logit.reg)
+lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)
-HLtest(object = fit)
+</code>
-# La valeur de p est de 0.9051814. Donc, on ne rejète pas notre modèle.
+++++ Réponse défi 5 |
-# L'ajustement du modèle est bon.
+<code rsplus | Coder MLMs>
-</file>
+m1 <- lmer(Z_TP~Z_Length + (1|Lake) + (1|Species), data = data, REML=TRUE)
+m2 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake) + (1+Z_Length|Species), data = data, REML=TRUE)
+m3 <- lmer(Z_TP~Z_Length + (1|Species), data = data, REML=TRUE)
+m4 <- lmer(Z_TP~Z_Length + (1|Lake), data = data, REML=TRUE)
+m5 <- lmer(Z_TP~Z_Length + (1+Z_Length|Species), data = data, REML=TRUE)
+m6 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake), data = data, REML=TRUE)
+m7 <- lmer(Z_TP~Z_Length + (1+Z_Length|Lake) + (1|Species), data = data, REML=TRUE)
+m8 <- lmer(Z_TP~Z_Length + (1|Lake) + (1+Z_Length|Species), data = data, REML=TRUE)
-** DÉFI 6 **
+# Modèle bonus!
+M0<-lm(Z_TP~Z_Length,data=data)
-Évaluez l'ajustement et le pouvoir prédictif du modèle ''model.bact2''.
+# Il est toujours utile de construire le modèle linéaire de base sans facteurs avec variation de l'ordonnée
-Comment pouvez-vous améliorer le pouvoir prédictif du modèle ?
+# à l'origine ou de pente pour voir la variation dans les valeurs de AICc values (même si "lm" n'utilise pas
+# la même méthode d'estimation).
-++++ Défi 6 : Solution |
+</code>
-<file rsplus | >
-null.d <- model.bact2$null.deviance
-resid.d <- model.bact2$deviance
-bact.pseudoR2 <- (null.d - resid.d) / null.d
-bact.pseudoR2
-# 0.0624257
-# C'est très faible !
-library(binomTools)
-HLtest(Rsq(model.bact2))
-# Chi-square statistic:  7.812347  with  8  df
-# P-value:  0.4520122
-# L'ajustement est adéquat.
-</file>
-Le pouvoir prédictif pourrait être augmenté en incluant plus de variables explicatives.
 ++++
-==== Représentation graphique des résultats ====
+----
-Lorsque le modèle a été validé, il peut être utile de représenter les résultats graphiquement afin de voir comment la variable est influencée par les variables explicatives. Une façon de faire est de mettre à la fois les valeurs observées et prédites de la variable réponses en fonction d'une variable explicative sur un même graphique. Voici un exemple avec le paquet ''ggplot2''. Revoyez [[r_workshop4|l'atelier 4]] pour plus d'informations sur ce paquet.
+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é.
-<file rsplus | Représenter graphiquement les résultats d'une régression logistique>
+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.
-library(ggplot2)
-ggplot(mites, aes(x = WatrCont, y = pa)) + geom_point() +
-stat_smooth(method = "glm", family= "binomial", se = FALSE) + xlab("Water content") +
-ylab("Probability of presence") +
-ggtitle("Probability of presence of Galumna sp. as a function of water content")
-</file>
-{{::logistic_regression_plot.png?500|}}
+<code rsplus | Comparer les MLMs>
-==== Un exemple avec des données de proportions ====
+# 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)
-Les données de proportions sont plus près des variables binaires que vous ne l'imaginez ! Supposons que vous êtes sur le terrain pour estimer la prévalence d'une maladiechez le cerf de Virginie le long d'un gradient environnemental de disponibilité en ressources. Vous échantillonnez dix individus dans dix populations différentes pour estimer la prévalence de la maladie. Votre feuille de données comporte dix lignes avec l'information suivante : code de la population et nombre d'individus infectés dans la population. À première vue, on dirait des données d'occurrences (//i.e.// count data), mais ce n'est pas correct de le voir ainsi : vous savez combien d'individus vous avez échantillonnés dans chaque population. Vous avez donc des données de proportions : c'est le nombre d'individus infectés sur dix individus échantillonnés. Si vous vous rappelez ce que vous avez lu dans la section sur les distributions, vous devez choisir une distribution binomiale pour modéliser ces données. Illustrons ceci avec un exemple :
-<file rsplus | Un GLM avec des données de proportions>
+# ii) Comparer les modèles en comparant les valeurs AICc
-# Commençons par générer des données en se basant sur l'exemple précédent :
+# Calculer les valeurs AICc pour chaque modèle
-# On choisit un nombre aléatoire entre 1 et 10 pour déterminer le nombre de cerfs infectés (objet 'n.infected').
+AICc<-c(AICc(M0), AICc(M1), AICc(M2), AICc(M3), AICc(M4), AICc(M5), AICc(M6), AICc(M7), AICc(M8))
-# On échantillonne dix individus dans dix populations différentes (objet 'n.total').
+# Mettre des valeurs dans une table pour faciliter la comparaison
-# 'res.avail' est un indice de qualité de l'habitat (disponibilité de ressources).
+Model<-c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8")
-set.seed(123)
+AICtable<-data.frame(Model=Model, AICc=AICc)
-n.infected <- sample(x = 1:10, size = 10, replace = TRUE)
+AICtable
-n.total <- rep(x = 10, times = 10)
+# M8 a la plus faible valeur AICc donc le meilleur modèle
-res.avail <- rnorm(n = 10, mean = 10, sd = 1)
+# M2 est également un bon modèle, mais tous les autres modèles ne sont pas aussi bons.
-# Ensuite, on spécifie le modèle. Prenez note comment on spécifie la variable réponse.
+</code>
-# Il faut indiquer le nombre de fois où la maladie a été détectée
-# et le nombre de fois où celle-ci n'a pas été détectée.
-prop.reg <- glm(cbind(n.infected, n.total - n.infected) ~ res.avail, family = binomial)
-summary(prop.reg)
-# Si vos données sont déjà sous la forme de proportions, voici comment faire dans R :
-# Créons tout d'abord un vecteur de proportions :
-prop.infected <- n.infected / n.total
-# On doit spécifier l'argument "weights" de la fonction glm() pour indiquer le nombre d'essais par site.
-prop.reg2 <- glm(prop.infected ~ res.avail, family = binomial, weights = n.total)
-summary(prop.reg2)
-# Les sommaires des modèles prop.reg et prop.reg2 sont identiques !
-</file>
-===== GLM avec des données d'abondance =====
+++++ Sortie|
+{{:fig_23_w5.png?200|}}
+++++
-Afin d'illustrer l'utilisation des GLMs avec des données d'abondance, nous allons utiliser un nouveau jeux de données; faramea.
+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.
-<file rsplus | loading data>
+<code rsplus | Recoder les MLMs>
-faramea <- read.csv(‘faramea.csv’, header = TRUE)
+# Une fois que les meilleurs modèles sont sélectionnés il faut remettre REML=TRUE
-</file>
+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)
+</code>
-Ce jeux de données s'intéresse à l'espèce d'arbre //Faramea occidentalis// sur l'île Barro Colorado au Panama. 43 transects ont été utilisés afin de mesurer le nombre d'arbre le long d'un gradient environnemental. Des caractéristiques environnementales, comme l'élévation ou la précipitation, ont aussi été mesurées au niveau de chaque transect. Examinons maintenant à quoi ressemble la distribution du nombre d'arbres par transect.
+----
+**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?
-<file rsplus | plot the data>
+++++ Réponse défi 6 |
-hist(faramea$Faramea.occidentalis, breaks=seq(0,45,1), xlab=expression(paste("Number of ",
+M2: La position trophique est une fonction de la longueur. L'intercept et l'effet de la longueur sur la position trophique peuvent varier selon l'espèce de poissons et le lac.
-italic(Faramea~occidentalis))), ylab="Frequency", main="", col="grey")
-</file>
+M8: La position trophique est une fonction de la longueur. L'intercept et l'effet de la longueur sur la position trophique peut varier selon l’espèce de poissons, mais seulement l'intercept peut varier par lac (et non la pente de la position trophique sur la longueur).
-{{ :plot_faramea_abd.png?500 |}}
+Biologiquement parlant, M2 indique que les facteurs intrinsèques des espèces (par exemple les taux de croissance) et des lacs  (par exemple, la productivité, la composition de la communauté, etc.) sont à la base de relations différentes entre la position trophique et la longueur (c.-à-d. pentes et intercepts), tandis que M8 indique que les facteurs intrinsèques des espèces seulement sont responsables pour les différentes relations (c.-à-d. pentes) et que, en moyenne, les positions trophiques pourraient être supérieures ou inférieures d’un lac à l’autre (par exemple intercepts).
-Nous pouvons remarquer qu'il n'y a que des valeurs entières et positives. Etant donné cette spécificité propore aux données de dénombrement, la distribution de Poisson, décrite ci dessus, semble un choix approprié pour modéliser ces données.
+Ces modèles sont très similaires dans leur structure et les unités AIC le suggèrent. La complexité supplémentaire de permettre que les pentes varient par lac dans M2 n’améliore pas le pouvoir prédictif par rapport au modèle M8.
-==== GLM avec une distribution de Poisson ====
-Comme nous l'avons vu plus haut, la distribution de Poisson est particulièrement appropriée pour modéliser des données de dénombrement car :
-  * elle ne spécifie des probabilités que pour des valeurs entières
-  * P(y<0) = 0, en d'autres termes la probabilité d'observer une valeur négative est nulle
-  * la relation entre la moyenne et la variance permet de manipuler des données hétérogènes (e.g. quand la variance dans les données augmente avec la moyenne)
-Dans notre exemple, nous voulons tester si l'élévation (une variable explicative continue) influence l'abondance de //Faramea occidentalis//. Il est généralement bien de commencer par utiliser un GLM et une distribution de Poisson lorsque nous cherchons à modéliser des données d'abondance. Nous allons donc commencer par ajuster ce type de modèle pour modéliser l'abondance de //Faramea occidentalis// en fonction de l'élévation.
+++++
-{{ ::plot_faramea_vs_elevation.png?400 |}}
+----
+====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.
+<code rsplus | 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)
+</code>
+++++ Sortie|
+{{:fig_24_w5.png?400|}}
+++++
-** Mais qu'est ce qui se cache derrière un tel GLM ?**\\
+L'étendue similaire des résidus suggère que le modèle est adéquat pour bien modéliser nos données.
-un GLM avec une distribution de Poisson considère que les valeurs y<sub>i</sub> de la variable réponse ont été générées par une distribution de Poisson dont la moyenne et la variance sont égales à //µ//<sub>i</sub>.
-Y<sub>i</sub> ∼ Poisson(//µ//<sub>i</sub>)
+Pour vérifier la supposition d'indépendance, il faut faire un graphique des résidus en fonction de chaque covariable du modèle :
-Avec E(Y<sub>i</sub>) = Var(Y<sub>i</sub>) = //µ//<sub>i</sub>
+<code rsplus | 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.
-Souvenez vous qu'un GLM est composé de deux parties, le prédicteur linéaire et la fonction de lien. Le prédicteur linéaire est utilisé comme une combinaison linéaire des différentes variables explicatives dont les effets sont estimés à travers les paramètres //β//. Le prédicteur linéaire peut être défini comme :
+# Espèce
+boxplot(E1 ~ Fish_Species,
+        ylab = "Normalized residuals",
+        data = data, xlab = "Species")
+abline(h = 0, lty = 2)
-//β//<sub>0</sub> + **X**<sub>i</sub>.//β//
+# Lac
+boxplot(E1 ~ Lake,
+        ylab = "Normalized residuals",
+        data = data, xlab = "Lake")
+abline(h = 0, lty = 2)
+</code>
+++++ Sortie|
+----
+{{:fig_25_w5.png?400|}}
+----
+{{:fig_26_w5.png?400|}}
+----
+{{:fig_27_w5.png?400|}}
+----
+++++
+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.
-où **X** est la matrice des variables explicatives, //β//<sub>0</sub> l'ordonnée à l'origine et //β// le vecteur des paramètres estimés.
+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.
-La fonction de lien, entre le prédicteur linéaire et la moyenne de la distribution //µ//<sub>i</sub> est la fonction logarithmique (relation log-linéaire) :
+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é.
-  * log(//µ//<sub>i</sub>) = //β//<sub>0</sub> + **X**<sub>i</sub>.//β//
+<code rsplus | Validation  du modèle >
-ou
+#D. Vérifier la normalité : histogramme
-  * //µ//<sub>i</sub> = exp(//β//<sub>0</sub> + **X**<sub>i</sub>.//β//)
+hist(E1)
+</code>
+++++ Sortie|
+{{:fig_28_w5.png?400|}}
+++++
-La distribution de Poisson donne la probabilité qu'une valeur Y<sub>i</sub> soit observée étant donné une moyenne //µ//<sub>i</sub> = exp(//β//<sub>0</sub> + **X**<sub>i</sub>.//β//). Dans ce modèle, les paramètres inconnus sont inclus dans le vecteur des coefficients de régression //β// (plus l'ordonnée à l'origine //β//<sub>0</sub>) et seront estimés à l'aide du maximum de vraisemblance.
+====2.4 Interpréter les résultats et les visualiser graphiquement ====
-Pour ajuster un GLM avec une distribution de Poisson sous R, il suffit de spécifier //family = poisson// dans la fonction glm(). Par défaut la fonction de lien est la fonction logarithmique.
-<file rsplus | fit a Poisson GLM>
+Vous pouvez voir le résumé du modèle à l'aide de :
-glm.poisson = glm(Faramea.occidentalis~Elevation, data=faramea, family=poisson)
+<code rsplus | Interprétation des résultats>
-summary(glm.poisson)
+# Vérifiez le résumé du modèle
-</file>
+# 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)
+</code>
+++++ Sortie|
+{{:fig_29_w5.png?500|}}
+++++
-Le résumé est similaire à celui de la fonction '**lm**' (voir atelier 3) et donne les estimations des paramètres. Vous pouvez aussi récupérer les estimations des paramètres à l'aide des fonctions suivantes :
+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) :
-<file rsplus | parameter estimates>
+{{:fig_30_w5.png?500|}}
-# ordonnée à l'origine
-summary(glm.poisson)$coefficients[1,1]
-# coefficient de regression de l'élévation
-summary(glm.poisson)$coefficients[2,1]
-</file>
+{{:fig_31_w5.png?500|}}
+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.
-==== La validation du modèle et le problème de la surdispersion ====
+{{:fig_32_w5.png?500|}}
-Un aspect important du résumé se trouve dans les dernières lignes :
+----
-<file rsplus | Poisson GLM summary>
+**DÉFI 7**\\
-summary(glm.poisson)
+a) Quelle est la pente et l'intervalle de confiance de la variable Z_Length du le modèle M8?
-# Null deviance: 414.81  on 42  degrees of freedom
-# Residual deviance: 388.12  on 41  degrees of freedom
-</file>
-L'estimation du maximum de vraisemblance est utilisé afin d'estimer les paramètres. Nous avons déjà mentionné que la déviance est l'équivalent en maximum de vraisemblance des sommes des carrés dans un modèle linéaire. Ici vous pouvez considérer la déviance nulle et la déviance résiduelle comme les équivalents de la somme totale des
+b) Est-ce que la pente de Z_Length est significativement différente de 0?
-carrés et de la somme des carrés résiduelle. La déviance résiduelle correspond à deux fois la différence entre la log-vraisemblance de deux modèles : un modèle qui s'ajuste parfaitement aux données (i.e. un modèle saturé) et le modèle que nous voulons tester. Si notre modèle est correct, la distribution de la déviance résiduelle est estimée selon une distribution du χ² avec //n//-//p//-1 degrés de liberté (où //n// correspond au nombre d'observations et //p// correspond au nombre de variables explicatives). Une implication très importante en ce qui nous concerne est que la déviance résiduelle doit être égale au nombre de degrés de liberté résiduels. Dans notre example, la déviance résiduelle vaut 388.12, tandis que nous avons 41 (43-1-1) degrés de liberté. La déviance est 9.5 fois supérieure au nombre de dégrés de liberté. Le modèle peut alors être qualifié de **surdispersé**.
-**La surdispersion**
+++++ Réponse défi 7 |
-La surdispersion peut être évaluée à l'aide du paramètre de surdispersion φ qui se mesure donc de la facon suivante :
+a) Quelle est la pente et son intervalle de confiance de la variable Z_Length dans le modèle M8?
+  - Pente = 0.4223
-                        φ = déviance résiduelle / dégrés de liberté résiduels
-  * φ < 1 indique qu'il y a sousdispersion
+<file>
-  * φ = 1 indique que la dispersion est conforme aux attendus
+   limite supérieure de l’IC = 0.4223 + 0.09*1.96 = 0.5987
-  * φ > 1 indicates qu'il y a surdispersion
+   limite inférieure de l’IC = 0.4223 - 0.09*1.96 = 0.2459
-Mais pourquoi un GLM présente-il de la surdispersion ? En fait, des modèles GLM sont surdispersés quand la variance dans les données est encore plus grande que ce qu'autorise la distribution de Poisson. Par exemple, cela peut se produire lorsque les données contiennent de nombreux zeros ou beaucoup de très grosses valeurs. Si nous revenons sur la distribution de nos données (ci dessus) nous pouvons remarquer que ces deux problèmes sont présents et que la distribution de Poisson n'était peut être pas le choix idéal. La surdispersion peut aussi survenir lorsque des variables explicatives ou des termes d'intéractions sont absentes ou bien encore lorsque qu'il y a des problèmes de valeurs aberrantes.
-La distribution de Poisson peut tenir compte de l'hétérogénéité présente dans des données grace à la relation entre sa moyenne et sa variance. Toutefois dans certains cas la variance augmente bien plus rapidement par rapport à la moyenne si bien que la distribution de Poisson n'est plus appropriée. Pour nous convaincre une dernière fois d'abandonner la distribution de Poisson pour modéliser l'abondance de l'espèce //faramea//  nous pouvons rapidement calculer la moyenne et la variance dans notre jeux de données :
-<file rsplus | mean and variance of the data>
-mean(faramea$Faramea.occidentalis)
-var(faramea$Faramea.occidentalis)
 </file>
-Dans la pratique, les GLMs basés sur la distribution de Poisson sont très pratique pour décrire la moyenne //µ//<sub>i</sub> mais vont sous-estimer la variance dans les données dès qu'il y a de la surdispersion. Par conséquent, les tests qui découlent du modèle seront trop laxistes. Il y a deux moyens de traiter les problèmes de surdispersion que nous allons détailler ci-dessous :
+b) Est-ce que la pente de Z_Length est significativement différente de 0?
-  * corriger la surdispersion en utilisant un **GLM quasi-Poisson **
+  - Oui, car l'IC n'inclut pas 0.
-  * choisir une nouvelle distribution comme la **binomiale négative**
+++++
+----
-==== GLM avec une correction "quasi" Poisson ====
+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 ()".
-Le principe d'un GLM avec une correction "quasi" Poisson est très simple; le paramètre de surdispersion (φ) est ajouté dans l'équation qui spécifie la variance du modèle :
-E(Y<sub>i</sub>) = //µ//<sub>i</sub>
+Deux façons de visualiser ces données sont :
+  - Montrer le modèle au niveau du groupe  (toutes les données groupées)
+  - Montrer le modèle au niveau de l’espèce ou du lac
-Var(Y<sub>i</sub>) = φ.//µ//<sub>i</sub>
+. Pour montrer le modèle au niveau du groupe :\\
-Le prédicteur linéaire, ainsi que la fonction de lien (log) restent les mêmes. La seule différence est que φ va être estimé afin de corriger le modèle. Les estimations des paramètres seront eux aussi inchangés, mais leurs écarts-types seront multiplés par √φ. Ainsi, certains paramètres qui étaient marginalement significatifs peuvent ne plus le rester.
+Obtenir les paramètres d'intérêts \\
-Dans R, la famille '**quasipoisson**' peut être utilisée pour traiter ces problèmes de surdispersion (de la même manière la famille **quasibinomial**' peut être utilisée). L'estimation de φ sera donné dans le résumé du modèle GLM quasi Poisson. Nous pouvons ajuster ce modèle de deux manières différentes :
+{{:fig_33_w5.png?600|}}
-<file rsplus | fit a new quasi-Poisson GLM>
+et tracer les données avec le modèle superposée
-# Option 1, nous ajustons un nouveau modèle GLM quasi-Poisson
-glm.quasipoisson = glm(Faramea.occidentalis~Elevation, data=faramea, family=quasipoisson)
-# Option 2, nous actualisons le modèle précédent :
-glm.quasipoisson = update(glm.poisson,family=quasipoisson)
-# regardons le résumé
-summary(glm.quasipoisson)
-</file>
-En examinant le résumé du modèle nous pouvons voir que φ est estimé à 15.97. Nous avons donc eu raison de corriger le modèle afin de prendre en compte la surdispersion. Par contre, si nous regardons la significativité du coefficient de regression associé à l'élévation, nous remarquons qu'il n'est plus significatif. Cependant, 15.97 indique que la surdispersion est forte et en général un GLM quasi-Poisson est favorisé lorsque φ est compris entre 1 et 15. Dans un soucis pédagogique, nous allons considérer que nous ne sommes pas satisfait avec ce modèle et que nous voulons maintenant ajuster une distribution binomiale négative à nos données.
+<code rsplus | 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.
-Deux points sont importants à garder en tête lorsque vous utilisez un GLM quasi-Poisson afin de corriger la surdispersion :
+# La première étape est d'obtenir les coefficients du modèle afin de les ajouter aux figures
+coef(M8)
-  * ** Les GLMs quasi-Poisson n'ont pas d'AIC.** En effet, la vraisemblance d'un modèle GLM quasi-Poisson ne peut pas être spécifiée et s'appuie sur une procédure de pseudo-maximum de vraisemblance. Par conséquence les GLMs quasi-Poisson n'ont pas d'AIC, et ce critère ne peut pas être utilisé afin de comparer différents modèles. Toutefois des alternatives ont été developpées pour gérer cette situation (e.g. quasi-AIC).
+# 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")
-  * ** La surdispersion influence la comparaison de modèles.** En effet, la surdispersion influence la comparaison de deux modèles emboités et doit donc être prise en considération. Par exemple, considérons que nous voulons comparer le modèle GLM1, qui contient //p//<sub>1</sub> paramètres avec le modèle GLM2, qui contient //p//<sub>2</sub> paramètres. GLM1 est emboité dans GLM2 et //p//<sub>2</sub> > //p//<sub>1</sub>. La comparaison des deux modèles est basées sur le test du rapport des vraisemblances des deux modèles, D<sub>1</sub> et D<sub>2</sub> respectivement. Si la surdispersion est connue, les déviances doivent être corrigées de manière approprié selon D* = D/φ, et le test final sera basé sur le critère D<sub>1</sub>* - D*<sub>2</sub> qui est supposé être distributé selon une distribution du χ² avec //p//<sub>2</sub>-//p//<sub>1</sub> degrés de liberté lorsque le modèle GLM1 est correct.
+Species.coef <- as.data.frame(coef(M8)$Fish_Species)
-  * Mais dans certain cas φ n'est pas connu. Par exemple, lorsque vous spécifiez un GLM avec un distribution normale. Dans ce cas, φ peut être estimé //a posteriori// en utilisant la déviance résiduelle du plus gros modèle de telle sorte que le critière de comparaison devienne [(D<sub>1</sub>-D<sub>2</sub>)/(//p//<sub>2</sub>-//p//<sub>1</sub>)]/[D<sub>2</sub>/(//n//-//p//<sub>2</sub>)]. Ce critère est supposé suivre une distribution F avec //p//<sub>2</sub>-//p//<sub>1</sub> et n-//p//<sub>2</sub> degrés de liberté.
+colnames(Species.coef) <- c("Intercept","Slope")
-==== GLM avec une distribution binomiale négative ====
-Un GLM avec une distribution binomiale négative (BN) est utilisé lorsque la surdispersion est très forte. La distribution BN contient un paramètre supplémentaire, //k//, qui va être très utile pour gérer les problèmes de surdispersion. Avant de rentrer dans les détails sur R, voyons rapidement ce qui se cache derrière la distribution BN. En fait, la distribution BN est la combinaison de deux distributions; une distribution de Poisson et une distribution Gamma. La distribution BN définie la distribution d'une variable aléatoire discrète de la même manière qu'une distribution de Poisson mais autorise la variance à être différente de la moyenne. Le mélange entre la distribution de Poisson et la distribution Gamma peut se résumer à l'aide de deux paramètres, //µ// et //k// qui spécifie la distribution de la facon suivante :
-Y ~ BN(//µ//, //k//)
-E(Y) = //µ// et Var(Y) = //µ// + //µ//²///k//
-De cette manière nous pouvons voir comment cette distribution va gérer la surdispersion dans les modèles GLM. Le deuxième terme de la variance de la distribution BN va déterminer le degré de surdispersion. En effet, la surdispersion est indirectement déterminée par //k//, que représente le paramètre de dispersion. Si //k// est grand (par rapport à //μ//²), la deuxième partie de la variance, //µ//²///k// va s'approcher de 0, et la variance de Y sera //μ//. Dans ce cas la distribution BN converge vers la distribution de Poisson et vous pourriez tout aussi bien utiliser cette dernière. Par contre, plus //k// sera petit et plus la surdispersion sera grande. Comme avec toutes les autres distributions, un GLM avec une distribution BN se spécifie en trois étapes. Tout d'abord le modèle fait l'hypothèse que les Y<sub>i</sub> suivent une distribution BN de moyenne //μ//<sub>i</sub> et de paramètre //k//.
-Y<sub>i</sub> ∼ BN(//µ//<sub>i</sub>, //k//)
-E(Y<sub>i</sub>) = //µ//<sub>i</sub> et Var(Y<sub>i</sub>) = //µ//<sub>i</sub> + //µ//<sub>i</sub>²///k//
-Les deux dernières étapes définissent le prédicteur linéaire ainsi que la fonction de lien entre la moyenne des Yi et le prédicteur linéaire. La fonction de lien utilisée par les GLMs avec une distribution BN est le logarithme ce qui permet de s'assurer que les valeurs prédites soient toujours positives.
-  * log(//µ//<sub>i</sub>) = //β//<sub>0</sub> + **X**<sub>i</sub>.//β//
-ou
-  * //µ//<sub>i</sub> = exp(//β//<sub>0</sub> + **X**<sub>i</sub>.//β//)
-La distribution NB n'est pas dans la fonction glm(), il faut donc installer et charger la librairie MASS. Si vous ne vous souvenez plus si cette librairie R est déjà installée, pas de problème utilisez la fonction suivante :
-<file rsplus | 'Mass' est il installé?>
-ifelse(length(which(installed.packages() == "MASS")) == 0,
-      {print("MASS not installed. Installing... "); install.packages("MASS")},
-      print("MASS already installed"))
-</file>
-Sinon, si vous savez que cette librairie n'est pas installée, vous pouvez directement taper la commande :
-<file rsplus | Installer 'MASS'>
-install.packages("MASS")
-</file>
-Souvenez vous de charger la librairie
-<file rsplus | Charger 'MASS'>
-library("MASS")
-</file>
-Vous pouvez ajuster un GLM avec une distribution BN à l'aide de la fonction glm.nb() :
-<file rsplus | Ajuster un GLM avec un BN'>
-glm.negbin = glm.nb(Faramea.occidentalis~Elevation, data=faramea)
-summary(glm.negbin)
-</file>
-Le résumé du modèle et similaire à celui des autres GLMs (e.g. GLMs Poisson). Cependant vous avez maintenant un nouveau paramètre, theta, qui est le paramètre //k// de la variance de votre distribution. L'écart-type de ce paramètre est aussi fourni, mais attention à son interprétation car l'intervalle n'est pas symétrique.
-==== Représentation graphique du modèle final ====
-Le GLM avec une distribution BN semble être le meilleur modèle pour modéliser nos données. Nous voulons maintenant représenter la relation entre le nombre de //Faramea occidentalis// et l'élévation.
-<file rsplus | Représenter le GLM'>
-# représenter les données observées
-plot(faramea$Elevation, faramea$Faramea.occidentalis, xlab="Elevation (m)", ylab=expression(paste("Number of", "  ", italic(Faramea~occidentalis))), pch=16, col=rgb(4,139,154,150,maxColorValue=255))
-# récupérons les valeurs de l'ordonnée à l'origine et de beta à partir du résumé et ajoutons les dans l'équation du modèle
-curve(exp(summary(glm.negbin)$coefficients[1,1]+summary(glm.negbin)$coefficients[2,1]*x),from=range(faramea$Elevation)[1],to=range(faramea$Elevation)[2],add=T, lwd=2, col="orangered")
-# récupérons les écarts-types pour construire l'intervalle de confiance du modèle
-curve(exp(summary(glm.negbin)$coefficients[1,1]+1.96*summary(glm.negbin)$coefficients[1,2]+summary(glm.negbin)$coefficients[2,1]*x+1.96*summary(glm.negbin)$coefficients[2,2]),from=range(faramea$Elevation)[1],to=range(faramea$Elevation)[2],add=T,lty=2, col="orangered")
-curve(exp(summary(glm.negbin)$coefficients[1,1]-1.96*summary(glm.negbin)$coefficients[1,2]+summary(glm.negbin)$coefficients[2,1]*x-1.96*summary(glm.negbin)$coefficients[2,2]),from=range(faramea$Elevation)[1],to=range(faramea$Elevation)[2],add=T,lty=2, col="orangered")
-</file>
-{{ :nb_glm_fit_to_faramea.png?450 |}}
-Nous pouvons voir que le nombre de //Faramea occidentalis// diminue de manière significative avec l'élévation. Toutefois, l'intervalle de confiance autour de notre modèle est assez large, notamment pour à faible élévation.
-==== Conclusion sur les GLM avec des données d'abondance ====
-Tous les GLMs que nous venons de voir pour modéliser des données de d'abondance (Poisson, quasi-Poisson et BN) utilisent la même relation log-linéaire entre moyenne et prédicteur linéaire (log(//µ//) = **X**.//β//). Toutefois ils vont autoriser différentes relations entre la moyenne et la variance et vont aussi se reposer sur des méthodes d'estimation de la vraisemblance différentes. Les GLMs Quasi-Poisson ou BN sont privilégiés afin de traiter la surdispersion. Malheureusement dans certains cas les données peuvent contenir trop de zeros et d'autres modèles seront plus éfficaces pour traiter ces situations. C'est par exemple le cas des //"zero-augmented models"// (e.g. zero-inflated Poisson [ZIP]) qui vont traiter les zéros indépendamment des autres valeurs.
-===== Les modèles linéaires généralisés mixtes (GLMM) =====
-On vient de voir que les variables réponses binaires et d'occurrences 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 [[r_atelier5|l'atelier 5]], 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 5) 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 [[http://www.facecouncil.org/puf/wp-content/uploads/2009/10/Generalized-linear-mixed-models.pdf|Bolker et al. (2009)]] et par la suite [[http://www.cell.com/cms/attachment/601623/4742452/mmc1.pdf|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).
-<file rsplus | 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")
-</file>
-=== 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.
-<file rsplus| 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)
-</file>
-{{ ::workshp_6_glmm_popu.png?600 |}}
-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.
-<file rsplus | 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")
-</file>
-{{ :workshp_6_glmm_boxplot.png?750 |}}
-{{ ::workshop_6_glmm_boxplot_popu.png?750 |}}
-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).
-  - Si nous observons une relation linéaire entre la variance et la moyenne avec une pente = 1, une famille de **Poisson** serait appropriée,
-  - 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,
-  - 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.
-{{ :worksp_6_error_dist.png?450 |}}
-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.
-<file rsplus | Poisson GLMM>
-mp1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
-               (1|popu)+
-               (1|gen),
-             data=dat.tf, family="poisson")
-</file>
-Nous pouvons vérifier la surdispersion en utilisant la fonction ''overdisp_fun'' fournie par [[http://www.cell.com/cms/attachment/601623/4742452/mmc1.pdf|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).
-<file rsplus | 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
-</file>
-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 **.
-<file rsplus | 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)
-</file>
-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 [[https://peerj.com/articles/616.pdf|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 ε<sub>i</sub> avec une distribution a priori lognormal. Une distribution Poisson-lognormal avec une moyenne µ et une variance a priori lognormal σ<sup>2</sup> est donné par :
-var(y) = µ + µ<sup>2</sup> [exp(σ<sup>2</sup>) - 1]
-En revanche, nous avons vu que la distribution binomiale négative (Poisson-gamma) était donné par:
-var(y) = µ + µ<sup>2</sup>/k
-En générale, la variance σ<sup>2</sup> 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.
-<file rsplus | 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)
-</file>
-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:
-  - l'approche théorique d'information (tel que le Critère d'Information d'Akaike; AIC), qui, comme nous l'avons vu dans [[r_workshop5#coding_potential_models_and_model_selection|l'atelier 5]] 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.
-  - 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).
-<file rsplus| É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
-</file>
-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:
-<file rsplus| 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")
-</file>
-{{::glmm_coefplot_int.png?450|}} {{::glmm_coefplot_pars.png?400|}}
-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)).
-<file rsplus| 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)
-</file>
-{{ ::glmm_dtoplot.png?450 |}}
-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|
-**Notez**: Pour cet ensemble de données, nous ne pouvons pas tester les pentes aléatoires, car lorsque les pentes aléatoires sont inclus dans le modèle, on obtient une forte corrélation entre les intercepts et pentes aléatoires. Ces fortes corrélations peuvent être dues à une "sur-paramétrisation" de notre modèle. Donc, bien qu'il pourrait y avoir une variation de l’effet nutriments ou herbivorie au niveau du génotype ou de la population (et bien que cette variation peut être ce qui nous intéresse le plus), il n’y a tout simplement pas suffisamment de données pour tester ces effets avec confiance.
-Regardons un exemple avec une pente aléatoire pour //clipping//:
-<file rsplus | Poisson-lognormal avec pentes aléatoires >
-# Poisson-lognormal avec pentes aléatoires pour popu et amd
-mpl2 <- glmer(total.fruits ~ nutrient*amd + rack + status +
-                (1|X) +
-                (amd|popu) +
-                (amd|gen),
-              data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))
-</file>
-Examinons les coefficients de variance en utilisant ''summary'' du modèle. Nous pourrions aussi examiner la matrice de corrélation entre les intercepts et pentes aléatoires. Une troisième option est d'utiliser la fonction ''printvc'' qui montre les matrices de variance-covariance avec un signe égal ( = ) à côté de toute composante de covariance avec une forte corrélation.
-<file rsplus | VarCorr>
-summary(mpl2) # option 1
-attr(VarCorr(mpl2)$gen,"correlation") # option 2
-printvc(mpl2) # option 3
-</file>
-Notez la corrélation parfaite entre (Intercept) et amdclipped (pente).
+# 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)
+</code>
+++++ Sortie|
+{{:fig_34_w5.png?500|}}
 ++++
+. Montrer les modèle au niveau de l’espèce ou du lac:\\
-=== Modèle final ===
+Obtenir les paramètres d'intérêts\\
-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.
+{{:fig_35_w5.png?200|}}
-<file rsplus| Évaluation des effets fixes avec AICc >
+et tracer les données avec le modèle superposé
-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"))
-</file>
-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).
+<code rsplus | 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")
-<file rsplus| Évaluation des effets fixes avec drop1>
+# Graphique 3 – Par Lac
-(dd_LRT <- drop1(mpl1,test="Chisq"))
+# Colorez les données par lac
-# Model:
+Plot_ByLake<-plot + geom_point(aes(colour = factor(Lake)), size = 4) + xlab("Length (mm)") + ylab("Trophic Position") + labs(title="By Lake") + fig
-#   total.fruits ~ nutrient * amd + rack + status + (1 | X) + (1 | popu) + (1 | gen)
+# Ajouter les lignes de régression avec les intercepts spécifiques à chaque lac
-#              Df    AIC   LRT   Pr(Chi)
+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")
-# <none>          5015.4
+</code>
-# rack          1 5070.5 57.083 4.179e-14 ***
+++++ Sortie|
-# status        2 5017.0  5.612   0.06044 .
+{{:fig_36_w5.png?500|}}
-# nutrient:amd  1 5016.8  3.444   0.06349 .
+{{:fig_37_w5.png?500|}}
-# ---
-#   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
-</file>
-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//).
-<file rsplus| 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)
-</file>
-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 7 **
-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:
-   - Quel est l'effet du type d'alimentation et du climat (effets fixes) sur PLD (variable réponse)?
-   - Est-ce que cette relation varie selon les taxons (effets aléatoires)?
-   - Quelle est la meilleure famille de distribution pour ces données?
-   - Une fois que vous avez déterminé la meilleure famille de distribution, réévaluez vos effets aléatoires et fixes.
-++++ Défi: Solution 1|
-<file rsplus | Effet du type d'alimentation et de la température sur PLD>
-# Charger les données
-inverts <- read.csv("inverts.csv")
-str(inverts)
-mlm1 <- lm(PLD ~ temp*feeding.type, data=inverts)
-anova(mlm1) # toutes les variables sont significatives
-</file>
 ++++
-++++ Défi: Solution 2|
+===== Conclusion et exercice de réflexion =====
-<file rsplus | Relations entre les taxons >
+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é.
-# Réponse vs effets fixes
+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.
-ggplot(inverts,aes(x=temp,y=log(PLD+1),colour=feeding.type)) +
-  geom_point() +
-  stat_summary(aes(x=as.numeric(temp)),fun.y=mean,geom="line") +
-  theme_bw() +
-  scale_color_manual(values=c("#3B9AB2","#F21A00")) + # from Wes Anderson Zissou palette
-  facet_wrap(~taxon)
-# Créer de nouvelles variables qui représentent toutes les combinaisons de feeding type x temp x taxa
-# (effets aléatoires)
-inverts <- within(inverts,
-{
-  # taxon x feeding.type
-  tft <- interaction(taxon,feeding.type,temp)
-  tft <- reorder(tft, PLD, mean)
-})
-# Boîte à moustaches total des fruits vs nouvelle variable (feeding type x temp x taxa)
+----
-ggplot(data = inverts, aes(factor(x = tft),y = log(PLD))) +
+**DÉFI 8**\\
-  geom_boxplot(colour = "skyblue2", outlier.shape = 21, outlier.colour = "skyblue2") +
+Situation :
-  theme_bw() + theme(axis.text.x=element_text(angle=90)) +
+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é.
-  stat_summary(fun.y=mean, geom="point", colour = "red")
-</file>
-++++
-++++ Défi: Solution 3|
+Quel modèle mixte pourriez-vous utiliser pour ce jeu de données?
-<file rsplus | Comparaison des familles de distribution >
-grpVars <- tapply(inverts$PLD, inverts$tft, var)
+++++ Réponse défi 8 |
-summary(grpVars)
+<code rsplus | Défi 8>
+>lmer(Bio_Div ~ Productivity + (1|Forest/Site))
-grpMeans <- tapply(inverts$PLD,inverts$tft, mean)
+># Ici, les effets aléatoires sont nichés (i.e. sites dans forêt) et non croisés.
-summary(grpMeans)
+</code>
-# Quasi-Poisson
-lm1 <- lm(grpVars~grpMeans-1)
-phi.fit <- coef(lm1)
-# Le -1 spécifie un modèle avec l'intercept fixer à zéro
-# Binomial négative
-lm2 <- lm(grpVars ~ I(grpMeans^2) + offset(grpMeans)-1)
-k.fit <- 1/coef(lm2)
-# Offset () est utilisé pour fixer la moyenne de chaque groupe à 1
-# Ajustement Loess non-paramétrique
-Lfit <- loess(grpVars~grpMeans)
-plot(grpVars ~ grpMeans, xlab="group means", ylab="group variances" )
-abline(a=0,b=1, lty=2)
-text(60,200,"Poisson")
-curve(phi.fit*x, col=2,add=TRUE)
-# bquote()est utilisé pour remplacer des valeurs numériques dans les équations avec des symboles
-text(60,800, bquote(paste("QP: ",sigma^2==.(round(phi.fit,1))*mu)),col=2)
-curve(x*(1+x/k.fit),col=4,add=TRUE)
-text(60,1600,paste("NB: k=",round(k.fit,1),sep=""),col=4)
-mvec <- 0:120
-lines(mvec,predict(Lfit,mvec),col=5)
-text(50,1300,"loess",col=5)
-# Poisson GLMM
-mp1 <- glmer(PLD ~ temp*feeding.type +
-               (1|taxon),
-             data=inverts, family="poisson")
-overdisp_fun(mp1)
-# rapport est significativement > 1
-# NB GLMM
-mnb1 <- glmer.nb(PLD ~ temp*feeding.type +
-                   (1|taxon),
-                 data=inverts)
-overdisp_fun(mnb1)
-# Semble bon!
-</file>
 ++++
+----
+----
+**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.
-++++ Défi: Solution 4|
+Quel modèle mixte pourriez-vous utiliser pour ce jeu de données?
-<file rsplus | Re-évaluer les effets aléatoires et fixes >
+++++ Réponse défi 9 |
+<code rsplus | Défi 9>
-# Re-évaluer les intercepts aléatoires
+> lmer(Mercury~Length*Habitat_Type+ (1|Site)...)
-summary(mnb1)$varcor
+</code>
-mnb1.taxless <- glm.nb(PLD ~ temp*feeding.type, data=inverts)
-# Ici, parce que nous comparons un glmer avec un glm, nous devons faire quelque chose de différent que
-# anova(). Pour tester l'importance de l’intercept aléatoire, nous allons comparer la vraisemblance de
-# chaque modèle :
-NL1 <- -logLik(mnb1)
-NL0 <- -logLik(mnb1.taxless)
-devdiff <- 2*(NL0-NL1)
-dfdiff <- attr(NL1,"df")-attr(NL0,"df")
-pchisq(devdiff,dfdiff,lower.tail=FALSE)
-# Nous pourrions aussi comparer l'AIC du modèle avec (mnb1) et sans (mnb1.taxless) effets aléatoires avec
-# la fonction AICtab()
-AICtab(mnb1,mnb1.taxless)
-# Changement important du AIC si nous supprimons l’intercept aléatoire. Donc, ça vaut la peine de garder
-# cet effet.
-# Graphique diagnostic
-locscaleplot(mnb1)
-# Graphique des paramètres de variance
-coefplot2(mnb1,ptype="vcov",intercept=TRUE, main="Random effect variance")
-# Graphique des effets fixes
-coefplot2(mnb1,intercept=TRUE,main="Fixed effect coefficient")
-# Graphique des intercepts aléatoires
-pp <- list(layout.widths=list(left.padding=0, right.padding=0))
-r2 <- ranef(mnb1,condVar=TRUE)
-d2 <- dotplot(r2, par.settings=pp)
-grid.arrange(d2$taxon,nrow=1)
-# Évaluer pentes aléatoires
-mnb2 <- glmer.nb(PLD ~ temp*feeding.type +
-                   (PLD|taxon),
-                 data=inverts)
-# Examiner composant de variance-covariance
-summary(mnb2) # option 1
-attr(VarCorr(mnb2)$taxon,"correlation") # option 2
-printvc(mnb2) # option 3
-# Forte corrélation entre les effets aléatoires -> pas assez de puissance pour tester pentes aléatoires
-# Re-évaluer les effets fixes
-# Remarque : pour utiliser la fonction de drop1 nous devons spécifier le paramètre thêta et exécuter
-# le modèle NB avec glmer :
-theta.mnb1 <- theta.md(inverts$PLD, fitted(mnb1), dfr = df.residual(mnb1))
-mnb1 <- glmer(PLD ~ temp*feeding.type +
-                (1|taxon),
-              data=inverts, family=negative.binomial(theta=theta.mnb1))
-(dd_LRT <- drop1(mnb1,test="Chisq"))
-(dd_AIC <- dfun(drop1(mnb1)))
-# Lorsque l’interaction feeding.type x température est supprimée, dAIC change de plus de 2 unité → suggère
-# de garder l'interaction dans le modèle
-</file>
 ++++
-===== Ressources additionnelles =====
+----
-**Livres**
+=====Livres de référence=====
+Gelman, A., and Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models (Cambridge University Press).
-B. Bolker (2009) [[http://ms.mcmaster.ca/~bolker/emdbook/|Ecological Models and Data in R]]. Princeton 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).
-A. Zuur et al. (2009) [[http://link.springer.com/book/10.1007/978-0-387-87458-6|Mixed Effects Models and Extensions in Ecology with R]]. Springer.\\
-**Sites web**
-[[http://glmm.wikidot.com/|GLMM for ecologists]]: Un site web très complet et avec une section Questions-Réponses pertinente !