Differences

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

--- r_atelier7 [2016/09/05 16:24]
vincent_fugere [Atelier 7: Modèles linéaires généralisés (mixtes)]
+++ r_atelier7 [2021/10/13 23:51] (current)
lsherin
@@ Line 1: / Line 1: @@
+<WRAP group>
+<WRAP centeralign>
+<WRAP important>
+<wrap em> __AVIS IMPORTANT__ </wrap>
+<wrap em> Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé. </wrap>
+<wrap em> Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le [[https://r.qcbs.ca/fr/workshops/r-workshop-06/|site web de la série d'ateliers R du CSBQ]]. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés. </wrap>
+<wrap em> Merci de votre compréhension, </wrap>
+<wrap em> Vos coordonnateurs de la série d’ateliers R du CSBQ. </wrap>
+</WRAP>
+</WRAP>
+<WRAP clear></WRAP>
 ======= Ateliers R du CSBQ =======
@@ Line 5: / Line 22: @@
 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 7: Modèles linéaires généralisés (mixtes) ======
+//Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil//
+<wrap em>AVIS IMPORTANT: MISES À JOUR MAJEURES</wrap>
+**Mise à jour de mars 2021:** Ce wiki n'est plus activement développé ou mis à jour. Le matériel mis à jour pour la série d'ateliers R du CSBQ est maintenant hébergé sur [[https://github.com/QCBSRworkshops/workshop06|la page GitHub]] des ateliers R du CSBQ.
+Le matériel disponible inclut;
+  - La [[https://qcbsrworkshops.github.io/workshop06/pres-fr/workshop06-pres-fr.html|présentation Rmarkdown]] pour cet atelier;
+  - Un [[https://qcbsrworkshops.github.io/workshop06/book-fr/workshop06-script-fr.R|script R]] qui suit la présentation.
+  - [[https://qcbsrworkshops.github.io/workshop06/book-fr/index.html|Le matériel écrit]] qui accompagne la présentation en format bookdown.
+**Mise à jour de janvier 2021:**
+Si vous cherchez l'atelier qui couvre les modèles linéaires généralisés, veuillez vous référer à l'**[[r_atelier7|Atelier 6]]**.
+Si vous cherchez l'atelier qui couvre les modèles linéaires et généralisés linéaires mixtes, veuillez vous référer à l'**[[r_atelier6|Atelier 7]]**.
+====== Atelier 6: Modèles linéaires généralisés ======
 Développé par : Cédric Frenette Dussault, Vincent Fugère, Thomas Lamy, Zofia Taranu
-**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.
+**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.
+**Lien vers la nouvelle [[https://qcbsrworkshops.github.io/workshop06/pres-fr/workshop06-pres-fr.html|présentation Rmarkdown]]**
+//S'il vous plaît essayez-la et dites aux coordonnateurs des ateliers R ce que vous en pensez!//
-Lien vers la présentation Prezi associée : [[https://prezi.com/ycczxaot6tre/|Prezi]]
+Lien vers l'ancienne [[https://prezi.com/ycczxaot6tre/|présentation Prezi]]
 Téléchargez le script R et les données pour cet atelier :
@@ Line 17: / Line 56: @@
   - [[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:
+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 4 et 6. À cet effet, chargeons des données et appliquons quelques modèles linéaires:
 <file rsplus | charger les données>
@@ Line 242: / Line 279: @@
 <file rsplus | Illustration du concept de prédicteur linéaire>
-# Chargez le jeu de données CO2. C'est ce jeu de données que nous avons utilisé dans l'atelier 3 !
+# Chargez le jeu de données CO2 utilisé dans un atelier précédent
 data(CO2)
 head(CO2)
@@ Line 281: / Line 318: @@
 </file>
-** DÉFI 5 **
+** DÉFI 1 **
 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 |
+++++ Défi 1 : Solution |
 <file rsplus | >
 model.bact1 <- glm(y ~ trt * week, family = binomial('logit'), data = bacteria)
@@ Line 397: / Line 434: @@
 </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 évaluer l'ajustement (//i.e.// goodness-of-fit) d'une régression logistique, les graphiques de diagnostic ne sont pas très utiles (voir atelier 4). 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'' ou ''vcdExtra''.
 <file rsplus | Le test de Hosmer-Lemeshow>
@@ Line 406: / Line 443: @@
 </file>
-** DÉFI 6 **
+** DÉFI 2 **
 Évaluez l'ajustement et le pouvoir prédictif du modèle ''model.bact2''.
 Comment pouvez-vous améliorer le pouvoir prédictif du modèle ?
-++++ Défi 6 : Solution |
+++++ Défi 2 : Solution |
 <file rsplus | >
 null.d <- model.bact2$null.deviance
@@ Line 429: / Line 466: @@
 ==== 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.
+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''. Revoir [[r_workshop3|l'atelier 3]] pour plus d'informations sur ce paquet.
 <file rsplus | Représenter graphiquement les résultats d'une régression logistique>
@@ Line 524: / Line 561: @@
 </file>
-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 :
+Le résumé est similaire à celui de la fonction '**lm**' (voir atelier 4) et donne les estimations des paramètres. Vous pouvez aussi récupérer les estimations des paramètres à l'aide des fonctions suivantes :
 <file rsplus | parameter estimates>
@@ Line 673: / Line 710: @@
-===== Les modèles linéaires généralisés mixtes (GLMM) =====
-On vient de voir que les variables réponses binaires et d'abondance peuvent être modélisées de façon plus appropriée en laissant tomber la supposition que les résidus doivent suivre une distribution normale en spécifiant une distribution pour les résidus (par exemple Poisson ou binomiale négative). Toutefois, que devrions-nous faire s'il existe une structure nichée dans les données ou des corrélations entre les observations qui causent certaines observations à être plus semblables les unes aux autres que les observations échantillonnées dans différents sites ou points dans le temps ? Comme nous l'avons vu dans [[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).
-++++
-=== Modèle final ===
-Nous terminons cette section sur les GLMMs avec l'évaluation des effets fixes. Nous décrivons d'abord les modèles potentiels et les comparons en utilisant AICc.
-<file rsplus| Évaluation des effets fixes avec AICc >
-mpl2 <- update(mpl1, . ~ . - rack) # modèle sans rack
-mpl3 <- update(mpl1, . ~ . - status) # modèle sans status
-mpl4 <- update(mpl1, . ~ . - amd:nutrient) # modèle sans l’ interaction amd:nutrient
-ICtab(mpl1,mpl2,mpl3,mpl4, type = c("AICc"))
-</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).
-<file rsplus| Évaluation des effets fixes avec drop1>
-(dd_LRT <- drop1(mpl1,test="Chisq"))
-# Model:
-#   total.fruits ~ nutrient * amd + rack + status + (1 | X) + (1 | popu) + (1 | gen)
-#              Df    AIC   LRT   Pr(Chi)
-# <none>          5015.4
-# rack          1 5070.5 57.083 4.179e-14 ***
-# status        2 5017.0  5.612   0.06044 .
-# nutrient:amd  1 5016.8  3.444   0.06349 .
-# ---
-#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-(dd_AIC <- dfun(drop1(mpl1)))
-#               Df   dAIC
-# <none>           0.000
-# rack          1 55.083
-# status        2  1.612
-# nutrient:amd  1  1.444
-</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|
-<file rsplus | Relations entre les taxons >
-# Réponse vs effets fixes
-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))) +
-  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>
-++++
-++++ Défi: Solution 3|
-<file rsplus | Comparaison des familles de distribution >
-grpVars <- tapply(inverts$PLD, inverts$tft, var)
-summary(grpVars)
-grpMeans <- tapply(inverts$PLD,inverts$tft, mean)
-summary(grpMeans)
-# 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: Solution 4|
-<file rsplus | Re-évaluer les effets aléatoires et fixes >
-# Re-évaluer les intercepts aléatoires
-summary(mnb1)$varcor
-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**
-B. Bolker (2009) [[http://ms.mcmaster.ca/~bolker/emdbook/|Ecological Models and Data in R]]. Princeton University Press.\\
-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 !