__AVIS IMPORTANT__ Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé. Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le [[https://r.qcbs.ca/fr/workshops/r-workshop-04/|site web de la série d'ateliers R du CSBQ]]. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés. Merci de votre compréhension, Vos coordonnateurs de la série d’ateliers R du CSBQ. ======= Ateliers R du CSBQ ======= [[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. //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// ====== Atelier 4 : Modèles linéaires ====== Développé par : Catherine Baltazar, Bérenger Bourgeois, Zofia Taranu, Shaun Turney, Willian Vieira **Résumé :** Dans cet atelier, vous apprendrez comment effectuer des modèles linéaires fréquemment utilisés en écologie tels que la régression simple, l’analyse de variance (ANOVA), l’analyse de covariance (ANCOVA) et la régression multiple avec le logiciel R. Après avoir vérifié les postulats de ces modèles (visuellement et statistiquement) et transformé vos données si nécessaire, l’interprétation des résultats et leur représentation graphique n’auront plus de secrets pour vous! **Lien vers la nouvelle [[https://qcbsrworkshops.github.io/workshop04/workshop04-fr/workshop04-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 l'ancienne [[http://prezi.com/tbpa702q4uxr/|présentation Prezi]] Téléchargez les scripts R et les données pour cet atelier : - [[http://qcbs.ca/wiki/_media/lm_f.r|Script R]] - [[http://qcbs.ca/wiki/_media/birdsdiet.csv| Données birdsdiet]] - [[http://qcbs.ca/wiki/_media/dickcissel.csv| Données dickcissel]] - [[http://qcbs.ca/wiki/_media/lm_adv.r|Script avancé (facultatif)]] - [[http://qcbs.ca/wiki/_media/my.plotcorr.r|Fonction Plotcorr (facultatif)]] - [[http://qcbs.ca/wiki/_media/col.plotcorr.r|Fonction Col plotcorr (facultatif)]] ===== Objectifs d'apprentissage ===== - Régression linéaire simple - Test de t - Analyse de la variance (ANOVA) - ANOVA à deux critères de classification - ANOVA non équilibrée\\ //(section avancée et optionnelle)// - Analyse de la covariance (ANCOVA) - Régression linéaire multiple ===== 1. Aperçu ===== ====1.1 Définir la moyenne et la variation==== En science, on s'intéresse souvent à déterminer les relations entre des variables. Dans cet atelier, nous apprendrons comment utiliser des modèles linéaires, un ensemble de modèles qui quantifie les relations entre des variables. On commencera par définir deux concepts-clés pour comprendre les modèles linéaires: la moyenne et la variation. La **moyenne** est une mesure de la valeur moyenne d'une population. Supposons que nous avons une variable aléatoire x, comme la grandeur des personnes dans une salle, et que nous voulons représenter les patrons de cette variable. En premier lieu, on peut utiliser la moyenne (qui peut être calculée de plusieurs manières). Par contre, la moyenne ne peut pas représenter une population au complet. On peut aussi décrire une population à l'aide de mesures de variation. La **variation** est la dispersion (ou l'écart) des observations autour de la moyenne. Par exemple, il y a peu de variation si tout le monde dans une salle ont presque la même grandeur, et beaucoup de variation si plusieurs personnes sont de grandeurs très différentes. L'écart moyen, la variance, l'écart type, et le coefficient de variation (ou l'écart type relatif) sont tous des mesures de variation que nous définirons ci-dessous. On peut mesurer l'écart de chaque élément par rapport à la moyenne: Avec l'écart de tous les éléments, on peut calculer l'**écart moyen**: Afin de transformer tous les variables en valeurs positives sans utiliser de valeurs absolues, on peut mettre chaque valeur au carré. On obtient alors la **variance**: Par contre, en mettant toutes les valeurs au carré, on change les unités de nos variables. Si on reprend notre exemple des grandeurs de personnes dans la salle, la variance sera en m², une unité qui n'est plus pertinente pour notre question initiale. Pour transformer ces valeurs en unités appropriées, on peut calculer l'**écart type**: Finalement, pour exprimer le **coefficient de variation** (également nommé l'écart type relatif) en pourcentage, nous avons: ====1.2 Les modèles linéaires==== Dans les modèles linéaires, on utilise les concepts de moyenne et de variation pour décrire la relation entre deux variables. On dit "modèles linéaires", parce qu'ils décrivent la relation entre variables avec une droite: {y_i} = {β_0} + {β_1}{x_{i1}} + ... + {β_p}{x_{ip}} + {ε_i} où\\ {y_i} est la variable réponse,\\ {β_0} est l'ordonnée à l'origine de la droite de régression,\\ {β_1} est le coefficient de variation de la 1^ère variable explicative,\\ {β_p} est le coefficient de variation de la p^ème variable explicative,\\ {x_i1} est la variable explicative continue pour la première observation,\\ {x_ip} est la variable explicative continue pour la p^ième observation,\\\\ {ε_i} sont les résidus du modèle (i.e. la variance inexpliquée).\\ La **variable réponse** (y) est la variable que nous voulons expliquer, ou la variable dépendante. Il n'y a qu'une variable réponse. Les **variables explicatives** (x) sont des variables qui peuvent (potentiellement) expliquer la variable réponse. Celles-ci sont les variables indépendantes. On peut inclure une ou plusieurs variables explicatives. Par exemple, supposons que nous voulons expliquer la variation en grandeur des personnes dans cette salle. La grandeur est la variable réponse. Des variables explicatives peuvent être l'âge ou le sexe. Dans les modèles linéaires, les variables réponses doivent être continues, mais les variables explicatives peuvent être continues ou catégoriques. Une **variable continue** a une infinité de valeurs possibles. Une **variable catégorique** a un nombre limité de valeurs possibles. L'âge, la température, et la latitude sont des variables continues. Le sexe, le stade de développement, et le pays sont des variables catégoriques. Pour les variables explicatives continues, le modèle linéaire évalue s'il y a une corrélation significative entre la variable réponse et la ou les variables explicatives. Pour les variables explicatives catégoriques, le modèle linéaire évalue si la valeur moyenne de la variable réponse diffère significativement entre les différents niveaux (ou groupes) des variables explicatives. Ceci deviendra plus clair en explorant les types de modèles linéaires dans les sections suivantes. Dans presque tous les cas, les variables explicatives n'expliquent pas toute la variation dans la variable réponse. Par exemple, le sexe et l'âge ne ne seront pas assez pour prédire la grandeur de chaque personne parfaitement. La variation qui reste inexpliquée sont les **résidus**, ou l'erreur. Le but d'un modèle linéaire est de trouver la meilleure estimation des paramètres (les variables β), puis d'évaluer la qualité de l'ajustement ("//goodness of fit//") du modèle. Plusieurs méthodes ont été développées pour calculer l'intercept et les coefficient de modèles linéaires. Le choix approprié dépend du modèle. Le concept général de ces méthodes consiste de minimiser les résidus. Selon le type et le nombre de variables explicatives considérées, différents outils statistiques peuvent être utilisés pour évaluer ces relations entre variables. Le tableau ci-dessous liste les 5 types d'analyses statistiques abordées dans cet atelier: ^ Analyse statistique^ Type de variable\\ réponse Y^ Type de variable\\ explicative X^ Nombre de\\ variables explicatives ^ Nombre de\\ niveaux k^ ^ Régression linéaire simple | Continue | Continue| 1| | ^ Test de t | ::: | Catégorique| 1| 2| ^ ANOVA | ::: | Catégorique| 1 (ANOVA à un facteur), 2\\ (ANOVA à deux facteurs) ou plus| 3 ou plus| ^ ANCOVA | ::: | Continue\\ ET catégorique| 2 ou plus| 2 ou plus| ^ Régression multiple | ::: | Continue| 2 ou plus| | ==== 1.3 Conditions de base du modèle linéaire ==== Pour être valide, les modèles linéaires s'appuient sur 4 conditions de base. Si les 4 conditions ne sont pas respectées, les résultats du modèle ne peuvent pas être interprétés de façon valable. - Les résidus sont **indépendants** - Les résidus suivent une **distribution normale** - Les résidus ont une **moyenne de 0** - Les résidus sont **homoscédastiques** (i.e. leur variance est constante) Notez que ces 4 conditions concernent les résidus, et non les variables réponses ou explicatives. Les résidus doivent être indépendants, c'est-à-dire qu'il n'y a pas de structure manquante dans le modèle (comme une autocorrélation spatiale ou temporelle). Les résidus doivent aussi suivre une distribution normale avec une moyenne de 0, signifiant que la majorité des résidus ont une valeur proche de 0 (i.e. l'erreur est très petite) et que la distribution est symmétrique (i.e. la variable réponse est sous-estimée autant qu'elle est surestimée). Les residus doivent être homoscédastiques, c'est-à-dire que l'erreur ne change pas beaucoup quand les variables explicatives changent de valeur. Dans les section suivantes, nous ne répétons pas les conditions ci-dessus pour chaque modèle. **Prenez conscience, par contre, que ces conditions de base s'appliquent à tous les modèles linéaires, incluant tous ceux qui seront abordés ci-dessous.** ====1.4 Statistiques de tests et p-values==== Une fois que le modèle a été exécuté dans R, on reçoit une **sortie du modèle** composé de plusieurs éléments. Comprendre chacun de ces éléments et identifier les éléments pertinents de la sortie demande un peu de pratique. La sortie inclut une estimation des paramètres (les variables β), ainsi que des statistiques de tests. Le statistique de test dépend du modèle linéaire employé (t est le statistique de test pour la régression linéaire et le test de t, et F est le statistique de test pour l'ANOVA). Dans un modèle linéaire, l'**hypothèse nulle** est typiquement qu'il n'y a aucune relation entre deux variables continues, out qu'il n'y a pas de différence entre les niveaux d'une variable catégorique. **Plus la valeur absolue d'un statistique de test est grande, moins il est probable que l'hypothèse nulle soit valide**. La probabilité exacte se trouve dans la sortie du modèle, et s'appelle le **"p-value"**. On peut penser au p-value comme la probabilité que l'hypothèse nulle est valide, malgré que c'est une simplification. //(Plus précisément, le p-value est la probabilité que, étant donné que l'hypothèse nulle est valide, la magnitude du statistique de test serait plus grande ou égale à le statistique de test réellement observé.)// Par convention, si un p-value est **inférieur à 0.05 (5%)**, l'hypothèse nulle est rejetée. Cette valeur de seuil s'appelle **α** (alpha). Si on rejette l'hypothèse nulle, on dit alors que l'hypothèse contraire est soutenue: il y a une relation significative entre variables ou une différence significative entre groupes. **Notez qu'on ne "prouve" pas d'hypothèses**, mais qu'on les soutient ou on les rejette. ==== 1.5 Flux de travail ==== Dans cet atelier, nous explorerons plusieurs types de modèles linéaires. La création et l'interprétation de chaque modèle diffère dans les détails, mais les principes généraux et le flux de travail s'appliquent à tous les types. Pour chaque modèle, on suivra donc les étapes de travail suivantes: - Visualiser les données (ceci peut aussi se faire plus tard) - Créer un modèle - Tester les 4 conditions de base du modèle - Ajuster le modèle si les conditions de base ne sont pas respectées - Interpréter les résultats du modèle ===== 2. Régression linéaire simple ===== La régression linéaire simple est un type de modèle linéaire qui contient **seulement une variable explicative continue**. La régression détermine si les deux variables (1 explicative, et 1 réponse) sont significativement corrélés. Une régression linéaire simple concerne deux paramètres qui doivent être estimés: l'**ordonnée à l'origine** (β₀) et un **coefficient de corrélation** (β₁). La méthode des moindres carrés est la méthode la plus couramment utilisée, et est employée par défaut dans la fonction ''lm()'' dans R. La méthode des moindres carrés fait passer une droite de manière à minimiser la somme des distances verticales au carré entre la droite et les données observées : autrement dit, la méthode vise à minimiser les résidus. Cliquez ci-dessous pour plus de détails mathématiques. À l'aide de la méthode des moindres carrés, le coefficient de corrélation (β₁) et l'ordonnée à l'origine (β₀) peuvent être calculés de la façon suivante : β_{1}={sum{i}{}{(x_{i}y_{i})}-overline{x}overline{y}}/sum{i}{}{(x_{i}-overline{x})}^2 = {Cov(x,y)}/{Var(x)} β_{0}=overline{y}-β_{1}overline{x} ==== 2.1 Effectuer un modèle linéaire ==== En utilisant le jeu de données ''bird'', nous allons exécuter une première régression linéaire sur l'abondance maximale en fonction de la masse. Dans R, une régression linéaire est implémentée à l'aide de la fonction ''lm()'' de la librairie ''stats'' : ''lm (y ~ x)'' |__**Note**__ : Avant d'utiliser une nouvelle fonction dans R, vous devriez vous référer à sa page d'aide (''?nomdelafonction'') afin de comprendre comment utiliser la fonction ainsi que les paramètres par défaut.| # Chargez les paquets et le jeu de données bird library(e1071) library(MASS) setwd("~/Desktop/...") # N'oubliez pas de spécifier votre répertoire de travail (note: le vôtre sera différent de celui-ci) bird<-read.csv("birdsdiet.csv") # Visualisez le tableau de données : names(bird) str(bird) head(bird) summary(bird) plot(bird) Le jeu de données bird contient sept variables: ^ Nom de la variable^ Description^ Type^ ^ Family | Nom de la famille | Chaînes de caractères| ^ MaxAbund | Abondance la plus élevée observée\\ en Amérique du Nord|Continue/numérique| ^ AvgAbund | Abondance moyenne sur tous les sites\\ en Amérique du Nord|Continue/numérique| ^ Mass | Taille corporelle en grammes| Continue/numérique| ^ Diet | Type de nourriture consommée| Catégorique – 5 niveaux (Plant; PlantInsect;\\ Insect; InsectVert; Vertebrate)| ^ Passerine| Est-ce un passereau?| Binaire (0/1)| ^ Aquatic | Est-ce un oiseau qui vit principalement dans ou près de l'eau?| Binaire (0/1)| Notez que Family, Diet, Passerine, et Aquatic sont tous des variables catégoriques, malgré le fait qu'ils soient encodés de façons différentes (chaîne de caractères, catégorique, binaire). Nous sommes maintenant prêts à exécuter le modèle linéaire : lm1 <- lm(bird$MaxAbund ~ bird$Mass) # où Y ~ X signifie Y "en fonction de" X> ==== 2.2 Validation des conditions de base ==== opar <- par(mfrow=c(2,2)) # Permet de créer les graphiques dans un panneau 2 x 2 plot(lm1) par(opar) # Remet la fenêtre graphique à un seul panneau === Vérifier l'indépendance === Les modèles linéaires s'appliquent seulement aux données indépendantes. En d'autres mots, le //y_i// correspondant à une certaine valeur //x_i// ne doit pas être influencée par d'autres valeurs //x_i//. Si vos données contiennent une certaine forme de structure de dépendance, comme une corrélation spatiale ou temporelle, l'hypothèse de l'indépendance est invalide. Malheureusement, il n'existe pas un graphique de diagnostique simple pour vérifier l'indépendance. Au contraire, il faut considérer son jeu de données avec soin et prudence. Est-ce qu'il y a un structure présente dans vos données qui crée une dépendance entre observations? Si les données s'agissent d'une série temporelle (donc, les observations proviennent des mêmes sites à plusieurs reprises), ou que plusieurs observations proviennent d'un même organisme, l'hypothèse d'indépendance est invalide. Il faut donc choisir un autre type de modèle. \\ ===Vérifier l'homoscédasticité et la moyenne résiduelle de 0=== //**Graphique des résidus vs. valeurs prédites**// - Le premier graphique de diagnostic est créé avec la fonction ''plot(lm1)''. Ce graphique illustre la dispersion des résidus en fonction des valeurs prédites par le modèle de régression linéaire. Chaque point représente la distance entre la variable réponse et la réponse prédite par le modèle. * Si les résidus sont dispersés de façon **aléatoire autour de la ligne de 0**, la relation est linéaire et la moyenne des résidus est 0. * Si les résidus forment une **bande horizontale** approximative autour de la ligne de 0, la variance des résidus est homogène (donc, ils sont homoscédastiques). * Si les résidus sont organisés **en forme d'entonnoir**, les résidus ne sont pas homoscédastiques. {{:workshop_3_lm1_residuals_vs_fitted.png?300|}} \\ //**Graphique "Scale-location"**// - Le troisième graphique de diagnostique permet de vérifier si la dispersion des résidus augmente pour une valeur prédite donnée (i.e. si la dispersion des résidus est causée par la variable explicative). Si la dispersion augmente, la condition de base d'homoscédasticité n'est pas respectée. {{:workshop_3_lm1_scale-location.png?300|}} \\ === Vérifier la normalité des résidus === //**Diagramme quantile-quantile (Q-Q)**// - La normalité des résidus peut être évaluée à l'aide d'un diagramme quantile-quantile. Ce graphique compare la distribution de probabilité des résidus du modèle à une distribution de probabilité de données normales. Si les résidus standardisés se trouvent près d'une ligne 1:1, les résidus peuvent être considérés comme normalement distribués. {{:workshop_3_lm1_qq.png?300|}} Dans ce cas-ci, les points ne sont pas bien alignés sur la droite, ce qui suggère que les résidus ne sont pas distribués normalement. \\ === Influence des observations aberrantes=== //**Diagramme de résidus vs. influence**// - En plus de valider les hypothèses de bases ci-dessus, on s'intéresse aussi à déterminer si certaines observations ont une forte influence. Bien qu'on ne teste pas un test de condition de base, ceci peut influencer notre interprétation des données. Si une ou certaines observations sont aberrantes (dont, si elles ont des valeurs très différentes des autres), le modèle peut être mal ajusté en raison de leur influence exagérée sur la calculation du modèle. Si (et seulement si!) ces observations correspondent à des erreurs de mesure ou à des exceptions, elles peuvent être retirées du jeu de données.\\ {{:workshop_3_lm1_leverage.png?300|}} ==== 2.3 Normalisation des données ==== Dans l'exemple précédent, les résidus du modèle ne suivaient pas une distribution normale, alors la condition de base de normalité est invalide. On peut quand même utiliser un modèle linéaire si on réussit à normaliser les données, afin de respecter la condition de normalité. L'étape suivante est donc de normaliser les données à l'aide de transformations mathématiques. Souvent, si on normalise les variables explicatives et/ou réponses, les résidus suivent une distribution normale. En plus des diagramme QQ, on peut évaluer la normalité d'une variable en traçant un histogramme avec la fonction ''hist()'', et en vérifiant visuellement que la variable suit une distribution normale. Par exemple : # Graphique Y ~ X et la ligne de régression plot(bird$MaxAbund ~ bird$Mass, pch=19, col="coral", ylab="Maximum Abundance", xlab="Mass") abline(lm1, lwd=2) ?plot # Pour obtenir plus de détails sur les arguments de la fonction plot(). # Allez voir colours() pour une liste de couleurs. # Les données sont-elles distribuées normalement ? hist(bird$MaxAbund,col="coral", main="Untransformed data", xlab="Maximum Abundance") hist(bird$Mass, col="coral", main="Untransformed data", xlab="Mass") {{:lm1_yvsx.png?300|}} {{:maxabund_hist.png?300|}} {{:mass_hist.png?300|}} Une troisième façon d'évaluer la normalité des données est d'utiliser le test de Shapiro-Wilk (fonction ''shapiro.test()''), qui compare la distribution des données observées à une distribution normale. Les hypothèses nulle et contraire sont: H₀: les données observées sont distribuées normalement\\ H₁: les données observées ne sont pas distribuées normalement Les données observées peuvent être considérées comme normalement distribuées lorsque la valeur de p calculée par le test de Shapiro-Wilk est supérieure au seuil α (généralement 0.05). # Teste l'hypothèse nulle que l'échantillon provient d'une population distribuée normalement shapiro.test(bird$MaxAbund) shapiro.test(bird$Mass) # Si p < 0.05, la distribution n'est pas normale # Si p > 0.05, la distribution est normale On peut également évaluer l'asymétrie d'une distribution avec la fonction ''skewness()'' : skewness(bird$MaxAbund) skewness(bird$Mass) # Une valeur positive indique une asymétrie vers la gauche (i.e. left-skewed distribution) # tandis qu'une valeur négative indique une asymétrie vers la droite (i.e. right skewed distribution). Les histogrammes, le test de Shapiro-Wilk, et le coefficient d'asymétrie ("skewness") indiquent tous que les variables doivent être transformées afin de respecter la condition de normalité (e.g. une transformation log₁₀). ==== 2.4 Transformation des données ==== Lorsque la condition de normalité n'est pas respectée, les variables peuvent être transformées afin d'améliorer la normalité de leur distribution en respectant ces règles : ^ Type de distribution^ Transformation^ Fonction R^ ^ Asymétrie positive modérée | sqrt{x} | sqrt(x)| ^ Asymétrie positive importante | log_10{(x)}|log10(x)| ^ Asymétrie positive importante | log_10{(x+C)}|log10(x + C) où C est une constante\\ ajoutée à chaque valeur de x afin que\\ la plus petite valeur soit 1| ^ Asymétrie négative modérée | sqrt{(K-x)}| sqrt(K - x) où K est une constante\\ soustraite de chaque valeur de x afin que\\ la plus petite valeur soit 1| ^ Asymétrie négative importante | log_10{(K-x)}| log10(K - x)| Dans notre cas, une transformation logarithmique (log₁₀) devrait être utilisée et enregistrée dans le tableau de données ''bird''. Le modèle peut ainsi être exécuté, vérifié, et interprété. # Ajoutez les variables transformées au tableau bird$logMaxAbund <- log10(bird$MaxAbund) bird$logMass <- log10(bird$Mass) names(bird) # pour visualiser le tableau avec les nouvelles variables hist(bird$logMaxAbund,col="yellowgreen", main="Log transformed", xlab=expression("log"[10]*"(Maximum Abundance)")) hist(bird$logMass,col="yellowgreen", main="Log transformed", xlab=expression("log"[10]*"(Mass)")) shapiro.test(bird$logMaxAbund); skewness(bird$logMaxAbund) shapiro.test(bird$logMass); skewness(bird$logMass) # Refaire l'analyse avec les transformations appropriées lm2 <- lm(bird$logMaxAbund ~ bird$logMass) # Reste-il des problèmes avec ce modèle (hétéroscédasticité, non-indépendance, forte influence)? opar <- par(mfrow=c(2,2)) plot(lm2, pch=19, col="gray") par(opar) {{:hist_logmaxabund.png?300|}} {{:hist_logmass.png?300|}} {{:plot_lm2_.png?550|}} ==== 2.5 Sortie du modèle ==== Une fois que les hypothèses (ou conditions) de base ont été validées, les résultats du modèle peuvent être interprétés. On obtient nos résultats avec la fonction ''summary()''. # Examinons les coefficients du modèle ainsi que les valeurs de p summary(lm2) # Vous pouvez faire apparaître seulement les coefficients lm2$coef # Quoi d'autre ? str(summary(lm2)) summary(lm2)$coefficients # où Std. Error est l'erreur type de chaque coefficient summary(lm2)$r.squared # Coefficient de détermination summary(lm2)$adj.r.squared # Coefficient de détermination ajusté summary(lm2)$sigma # Erreur type résiduelle (racine du carré moyen de l'erreur) # etc. # Vous pouvez également vérifier l'équation du coefficient de détermination (R2)par vous-mêmes : SSE <- sum(resid(lm2)^2) SST <- sum((bird$logMaxAbund - mean(bird$logMaxAbund))^2) R2 <- 1 - ((SSE)/SST) R2 La sortie de cette fonction présente tous les résultats du modèle : lm(formula = logMaxAbund ~ logMass, data = bird) Residuals: Min 1Q Median 3Q Max -1.93562 -0.39982 0.05487 0.40625 1.61469 Estimate Std. Error t value Pr(>|t|) (Intercept) 1.6724 0.2472 6.767 1.17e-08 *** logMass -0.2361 0.1170 -2.019 0.0487 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6959 on 52 degrees of freedom Multiple R-squared: 0.07267, Adjusted R-squared: 0.05484 F-statistic: 4.075 on 1 and 52 DF, p-value: 0.04869 Les **coefficients de régression** du modèle et leur **erreur type** se trouvent dans les deuxième et troisième colonnes du tableau de la régression, respectivement. Donc,\\ β₀ = 1.6724 ± 0.2472 est l'ordonnée à l'origine (± e.t.) du modèle de régression,\\ β₁ = -0.2361 ± 0.1170 est la pente (± e.t.) du modèle de régression. et finalement : //logMaxAbund// = 1.6724 (± 0.2472) - 0.2361 (± 0.1170) x //logMass// représente le modèle paramétré. Les valeurs de **t** et leurs **p-values** (dans les 4^ème et 5^ème colonnes du tableau, respectivement) testent si les coefficients de régression diffèrent significativement de zéro. Dans notre cas, on voit que la variable //logMass// a une influence significative sur la variable //logMaxAbund// parce que le p-value associée à au coefficient de variation (i.e. la pente) du modèle de régression est inférieure à 0.05. De plus, la relation entre ces deux variables est négative, car le coefficient de variation du modèle a une valeur négative.\\ **Rappel:** Une corrélation entre deux variables n'implique pas nécessairement une relation de causalité. Inversement, l'absence de corrélation entre deux variables n'implique pas nécessairement une absence de relation entre ces deux variables; c'est le cas, par exemple, lorsque la relation n'est pas linéaire. L'ajustement d'un modèle de régression linéaire est donné par le **R² ajusté** (ici 0.05484). Cette valeur mesure la proportion de variation qui est expliquée par le modèle. Pour plus de détails à propos de son calcul, cliquez ci-dessous: overline{R}^2=1-(1-R^2){n-1}/{n-p-1} où\\ p est le nombre total de paramètres de régression et n est la taille d'échantillon,\\ R^2={SS_reg}/{SS_tot} \\ {SS_tot}=sum{i}{}{({y_i}-overline{y})}^2 est la dispersion totale,\\ {SS_reg}=sum{i}{}{(hat{y_i}-overline{y})}^2 est la dispersion de la régression - aussi appelée la variance expliquée par le modèle. \\ Le R² ajusté varie entre 0 et 1. Plus ce coefficient est élevé (donc, plus qu'il s'approche de 1), plus le modèle est bien ajusté aux données. Dans ce cas-ci, la relation entre les variables //logMaxAbund// et //logMass// est très faible.\\ La dernière ligne de la sortie du modèle représente la **statistique F** du modèle et le **p-value** qui y est associée. Si la valeur de p est inférieure à 0.05, le modèle de régression décrit mieux la relation entre les variables qu'un modèle nul.\\ ==== 2.6 Visualisation graphique ==== Les résultats d'une régression linéaire sont généralement illustrés par un graphique de la variable réponse en fonction des variables explicatives. Une droite de régression y est tracée (et, si nécessaire, les intervalles de confiance) avec le code R suivant : plot(logMaxAbund ~ logMass, data=bird, pch=19, col="yellowgreen", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)")) abline(lm2, lwd=2) # On peut faire ressortir les points avec une forte influence points(bird$logMass[32], bird$logMaxAbund[32], pch=19, col="violet") points(bird$logMass[21], bird$logMaxAbund[21], pch=19, col="violet") points(bird$logMass[50], bird$logMaxAbund[50], pch=19, col="violet") # On peut également tracer les intervalles de confiance confit<-predict(lm2,interval="confidence") points(bird$logMass,confit[,2]) points(bird$logMass,confit[,3]) {{:yvxx_lm2.png?400|}} ==== 2.7 Sous-ensembles d'observations ==== Il est aussi possible d'analyser seulement une partie des observations. Par exemple, on peut refaire l'analyse de régression sur seulement les oiseaux terrestres. # Souvenez-vous qu'on peut exclure des valeurs avec le symbole "!" # On peut analyser un sous-ensemble des données de "bird" en utilisant l'argument 'subset' de la fonction lm(). lm3 <- lm(logMaxAbund ~ logMass, data=bird, subset =! bird$Aquatic) # enlever les oiseaux aquatiques du modèle # Cette commande permet également d'exclure les oiseaux aquatiques lm3 <- lm(logMaxAbund ~ logMass, data=bird, subset=bird$Aquatic == 0) # Examinons le modèle opar <- par(mfrow=c(2,2)) plot(lm3) summary(lm3) par(opar) # Comparons les deux analyses opar <- par(mfrow=c(1,2)) plot(logMaxAbund ~ logMass, data=bird, main="All birds", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)")) abline(lm2,lwd=2) plot(logMaxAbund ~ logMass, data=bird, subset=!bird$Aquatic, main="Terrestrial birds", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"), pch=19) abline(lm3,lwd=2) opar(par) ---- ---**DÉFI 1**--- Examinez la relation entre //log₁₀(MaxAbund)// et //log₁₀(Mass)// pour les passereaux (i.e. passerine birds).\\ ''Conseil :'' La variable 'Passerine' est codée comme 0 et 1, comme la variable 'Aquatic'. Vous pouvez vérifier ceci avec la commande ''str(bird)''. ++++ Défi 1 : Solution | lm4 <- lm(logMaxAbund ~ logMass, data=bird, subset=bird$Passerine == 1) # Examinez les graphiques de diagnostic opar <- par(mfrow=c(2,2)) plot(lm4) summary(lm4) par(opar) # Comparez la variance expliquée par lm2, lm3 et lm4 str(summary(lm4)) # Rappelez-vous qu'on veut le R^2 ajusté summary(lm4)$adj.r.squared # R2-adj = -0.02 summary(lm2)$adj.r.squared # R2-adj = 0.05 summary(lm3)$adj.r.squared # R2-adj = 0.25 # Comparez visuellement les trois modèles opar <- par(mfrow=c(1,3)) plot(logMaxAbund ~ logMass, data=bird, main="All birds", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"), pch=19, col="yellowgreen") abline(lm2,lwd=2) plot(logMaxAbund ~ logMass, subset=Passerine == 1, data=bird, main="Passerine birds", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"), pch=19, col="violet") abline(lm4,lwd=2) plot(logMaxAbund ~ logMass, data=bird, subset=!bird$Aquatic, main="Terrestrial birds", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"), pch=19, col="skyblue") abline(lm3,lwd=2) par(opar) **Conclusion :** Le meilleur modèle parmi les trois est lm3 (celui effectué sur seulement les oiseaux terrestres ) {{:lm2_lm3_lm4.png?850|}} ++++ ---- ===== 3. ANOVA ===== L'analyse de la variance (ANOVA) est un type de modèle linéaire pour **une variable réponse continue**, et **une ou plusieurs variables explicatives catégoriques**. Les variables explicatives catégoriques peuvent comprendre plusieurs niveaux (ou groupes). Par exemple, une variable "couleur" peut avoir 3 niveaux: vert, bleu, et jaune. L'ANOVA teste si la moyenne de la variable réponse diffère entre ces niveaux ou groupes. Par exemple, si la masse des bleuets diffère selon leur couleur. Le calcul de l'ANOVA se base sur la partition de la somme des carrés, et compare la **variance entre les traitements** à celle à l'intérieur des traitements (i.e. la **variance intra-traitement**). Si la variance entre les traitements est supérieure à la variance intra-traitement, la variable explicative a un effet plus important que l'erreur aléatoire (due à la variance intra-traitement). La variable explicative est donc susceptible d'influencer significativement la variable réponse. Dans un ANOVA, la comparaison de la variance entre les traitements à celle intra-traitement se fait en calculant la **statistique F**. Cette statistique correspond au ratio entre la moyenne des carrés des traitements (MS_Trt) et la moyenne des carrés des erreurs (MS_E). Ces deux termes sont obtenus en divisant leurs sommes des carrés respectives par leurs degrés de liberté, comme on le voit présenté typiquement dans un tableau d'ANOVA (cliquez pour voir un tel tableau ci-dessous). Un **p-value** peut ensuite être calculée à partir de la statistique de F, qui suit une distribution de khi carré (χ²) Cliquez pour plus de détails mathématiques. ^ Source de\\ variation^ Degrés de\\ liberté (df)^ Somme des carrés\\ des écarts à la moyenne^ Moyenne des carrés ^ Statistique de F^ ^ Total |ra-1| {SS_Tot}=sum{i,j}{}({y_ij}-overline{y})^2 | | | ^ Facteur A|a-1| {SS_FacteurA}= r sum{i}{}({overline{y}_i}-overline{y})^2|{MS_FacteurA}={SS_FacteurA}/{a-1}|F={MS_FacteurA}/{MS_E}| ^ Résidus |a(r-1)| {SS_E}= sum{i,j}{}({y_ij}-{overline{y}_i})^2 |{MS_E}={SS_E}/{a(r-1)}| | a: nombre de niveaux de la variable explicative A; r: nombre de répétitions par traitement; overline{y}: moyenne globale de la variable réponse; overline{y}_i : moyenne de la variable réponse du traitement i. ==== 3.1 Types d'ANOVA ==== - **ANOVA à un critère de classification**\\ Une variable explicative catégorique avec au moins 2 niveaux. S'il n'y a que 2 niveaux, un test de t peut être utilisé. - **ANOVA à deux critères de classification** //[[http://qcbs.ca/wiki/r_atelier4#anova_a_deux_criteres_de_classification|(voir la section ci-dessous)]]//\\ - Deux variables explicatives catégoriques ou plus,\\ - Chaque facteur peut avoir plusieurs niveaux,\\ - Les interactions entre chaque variable explicative catégorique doivent être testées. - **Mesures répétées**\\ L'ANOVA peut être utilisée pour des mesures répétées, mais ce sujet n'est pas abordé dans cet atelier. Un **modèle linéaire mixte** peut également être utilisé pour ce type de données //[[http://qcbs.ca/wiki/r_atelier6|(voir l'atelier 6)]]//. ==== 3.2 Test de t ==== Si on a une variable explicative ayant seulement **2 niveaux**, on peut utiliser un test de t de Student pour déterminer s'il y a une différence entre la moyenne des 2 niveaux. Selon les données, on peut choisir de tester une **hypothèse unilatérale**: c'est-à-dire qu'on peut déterminer si la moyenne d'un groupe est plus grande que celle de l'autre groupe, plutôt que de seulement détecter une différence entre les moyennes du groupe. Cliquez ci-dessous pour plus de détails mathématiques. Pour le test de t, la statistique t qui est utilisée pour déterminer la valeur de p est calculée de la manière suivante :\\ t= (overline{y}_1-overline{y}_2)/sqrt{{s_1}^2/n_1 + {s_2}^2/n_2} où overline{y}₁ et overline{y}₂ sont les moyennes de la variable réponse y pour les groupes 1 et 2 respectivement,\\ s₁² et s₂² sont les variances de la variable réponse y pour les groupes 1 et 2 respectivement,\\ n₁ et n₂ sont les tailles d'échantillons des groupes 1 et 2 respectivement. Notez que le test de t est mathématiquement équivalent à une ANOVA à un critère de classification ayant deux niveaux. \\ === Conditions de base === Si les conditions de base du test de t ne sont pas respectées, les résultats du test peuvent être erronés. Voici quelques notes à propos de la validation de ces conditions de base: - **Normalité des données**\\ Comme pour la régression linéaire simple, la variable réponse doit être distribuée normalement. Si cette condition n'est pas respectée, mais que la distribution est relativement symétrique, que la moyenne est près du centre de la distribution, et que la distribution est unimodale, le test de t donnera un résultat valable en autant que la taille de l'échantillon soit suffisante (règle empirique : ~30 observations). Si les données sont fortement asymétriques, il est nécessaire d'avoir un très large échantillon pour que le test fonctionne. Dans ce cas-là, il est préférable d'utiliser un test non-paramétrique. - **Homoscédasticité**\\ Une autre supposition importante du test de t est que les variances des deux groupes sont égales. Ceci permet de calculer une variance combinée qui est utilisée pour calculer l'erreur type de la différence des moyennes. Si les variances des deux groupes sont inégales, la probabilité de commettre une erreur de type I (i.e. rejeter l'hypothèse nulle alors qu'elle est vraie) est supérieure au seuil α.\\ La robustesse du test de t augmente avec la taille de l'échantillon et est supérieure lorsque les groupes sont de même taille.\\ Il est possible d'évaluer la différence de variance entre deux échantillons en se demandant quelle est la probabilité de tirer deux échantillons d'une population avec des variances identiques alors que les échantillons ont des variances de s₁² et s₂². \\ Pour ce faire, il faut effectuer un test de ratio des variances (i.e. un test de F). === Non-respect des conditions de base === Si les variances entre les groupes ne sont pas égales, il est possible de corriger la situation avec la [[http://fr.wikipedia.org/wiki/Test_t_de_Welch|correction de Welch]]. Si les conditions de base ne sont toujours pas respectées, il faut utiliser la version non paramétrique du test de t : [[http://en.wikipedia.org/wiki/Mann–Whitney_U_test|le test de Mann-Whitney]]. Finalement, si les deux groupes ne sont pas indépendants (e.g. mesures prises sur un même individu à deux périodes différentes), il faut utiliser un [[http://en.wikipedia.org/wiki/Student's_t-test#Paired_samples|test de t apparié]]. === Effectuer un test de t === Dans R, le test de t est exécutés avec la fonction ''t.test''. Par exemple, pour évaluer la différence de masse entre des oiseaux aquatiques et non aquatiques, on peut utiliser le script suivant :\\ # Test de t boxplot(logMass ~ Aquatic, data=bird, ylab=expression("log"[10]*"(Bird Mass)"), names=c("Non-Aquatic","Aquatic"), col=c("yellowgreen","skyblue")) # Tout d'abord, vérifions si les variances de chaque groupe sont égales # Note : il n'est pas nécessaire de vérifier la normalité des données, # car on utilise déjà une transformation logarithmique tapply(bird$logMass,bird$Aquatic,var) var.test(logMass~Aquatic,data=bird) # Nous sommes prêts pour le test de t ttest1 <- t.test(logMass~Aquatic, var.equal=TRUE, data=bird) # Cette commande est équivalente : ttest1 <- t.test(x=bird$logMass[bird$Aquatic==0], y=bird$logMass[bird$Aquatic==1], var.equal=TRUE) ttest1 {{:lm_fig_10.png?400|}} Two Sample t-test data: logMass by Aquatic t = -7.7707, df = 52, p-value = 2.936e-10 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.6669697 -0.9827343 sample estimates: mean of x mean of y 1.583437 2.908289 Ici, on voit que le test de ratio des variances n'est pas statistiquement différent de 1, ce qui signifie que les variances entre les groupes sont égales. Étant donné que notre valeur de p est inférieure à 0.05, l'hypothèse nulle (i.e. l'absence de différence de masse entre les deux groupes) est rejetée. === Test de t unilatéral === L'argument //"alternative"// de la fonction ''t.test()'' permet d'effectuer un test de t unilatéral. Par exemple, si on veut tester l'hypothèse que les oiseaux terrestres sont plus légers que les oiseaux aquatiques, on peut écrire la commande de la façon suivante : # Test de t en spécifiant l'argument "alternative" uni.ttest1 <- t.test(logMass~Aquatic, var.equal=TRUE, data=bird, alternative="less") uni.ttest1 Dans la sortie R obtenue par ''uni.ttest1'', les résultats du test sont indiqués à la troisième ligne : Two Sample t-test data: logMass by Aquatic t = -7.7707, df = 52, p-value = 1.468e-10 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf -1.039331 sample estimates: mean in group 0 mean in group 1 1.583437 2.908289 Dans ce cas-ci, la statistique du test de t est t = -7.7707 avec df = 52 degrés de liberté, ce qui donne une valeur de p = 1.468e-10. On rejette donc l'hypothèse nulle. On peut conclure que les oiseaux aquatiques sont significativement plus lourds que les oiseaux terrestres. === Effectuer un test de t avec lm()=== Un test de t est un modèle linéaire et un cas spécifique d'ANOVA avec une variable explicative ayant dexu niveaux. On peut alors effectuer un test de t avec la fonction ''lm()'' dans R. ttest.lm1 <- lm(logMass ~ Aquatic, data=bird) anova(ttest.lm1) Analysis of Variance Table Response: logMass Df Sum Sq Mean Sq F value Pr(>F) Aquatic 1 19.015 19.0150 60.385 2.936e-10 *** Residuals 52 16.375 0.3149 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Quand les variances sont équivalents, on peut montrer que t² = F: ttest1$statistic^2 anova(ttest.lm1)$F ==== 3.3 Effectuer une ANOVA ==== Le test de t s'applique seulement quand on a une seule variable explicative catégorique, qui comprend 2 niveaux. Pour tous les autres modèles linéaires avec des variables explicatives catégoriques, on utilise une ANOVA. Commençons tout d'abord par visualiser les données avec la fonction ''boxplot()''. Rappelez-vous que, dans R, les groupes sont ordonnés par ordre alphabétique par défaut. Il est possible de réorganiser les groupes autrement. Par exemple, on peut les ordonner par ordre croissant de la médiane de chaque diète.\\ Une autre façon de visualiser les effets des facteurs est d'utiliser la fonction ''plot.design()''. Cette fonction permet de représenter les valeurs moyennes des niveaux d'un facteur (par une ligne verticale) et la moyenne globale de la variable réponse (par une ligne horizontale). # Ordre alphabétique par défaut boxplot(logMaxAbund ~ Diet, data=bird) # Réorganiser l'ordre des facteurs med <- sort(tapply(bird$logMaxAbund, bird$Diet, median)) boxplot(logMaxAbund ~ factor(Diet, levels=names(med)), data=bird, col=c("white","lightblue1", "skyblue1","skyblue3","skyblue4")) plot.design(logMaxAbund ~ Diet, data=bird, ylab = expression("log"[10]*"(Maximum Abundance)")) {{:lm_fig_11.png?450|}} {{:plot_design.png?350|}} Nous sommes maintenant prêts à effectuer une ANOVA. Dans R, la fonction ''aov()'' permet d'effectuer une ANOVA directement. Il est également possible d'effectuer une ANOVA avec la fonction ''anova()'' qui exécute l'ANOVA comme un modèle linéaire : # En utilisant aov() aov1 <- aov(logMaxAbund ~ Diet, data=bird) summary(aov1) # En utilisant lm() anov1 <- lm(logMaxAbund ~ Diet, data=bird) anova(anov1) ==== 3.4 Vérification des conditions de base ==== Comme la régression linéaire simple et le test de t, l'ANOVA doit respecter 4 conditions statistiques pour que les résultats soient valides. Voici quelques conseils pour tester ces conditions pour une ANOVA: - **Distribution normale**\\ Les résidus d'un modèle d'ANOVA peuvent être visualisés à l'aide d'un diagramme quantile-quantile (Q-Q). Les résidus sont considérés comme normalement distribués s'ils se répartissent le long de la droite 1:1. Si ce n'est pas le cas, les résultats de l'ANOVA ne peuvent pas être interprétés. - **Homoscédasticité**\\ L'ANOVA est valide seulement lorsque la variance des résidus est homogène entre les groupes. Cette homoscédasticité peut être vérifiée par un graphique des résidus en fonction des valeurs prédites ou par un diagramme diagnostic "scale-location". Si ces graphiques montrent une dispersion équivalente des résidus pour chaque valeur prédite, la variance des résidus peut être considérée homogène.\\ Il est également possible d'effectuer un test de Bartlett à l'aide de la fonction ''bartlett.test()''. Si la valeur de p de ce test est supérieure à 0.05, l'hypothèse nulle H₀: s₁² = s₂² =... = s_j² =... = s_n² est acceptée (i.e. l'homoscédasticité est respectée).\\ Une transformation de la variable réponse peut être utilisée si cette supposition n'est pas respectée. - **Additivité**\\ Les effets de deux facteurs sont additifs si l'effet d'un facteur demeure constant pour tous les niveaux d'un autre facteur. Chaque facteur doit influencer la variable réponse de manière indépendante des autres facteurs. Si les conditions ne sont pas respectées, vous pouvez essayer de transformer la variable réponse. Ceci peut aider à normaliser les résidus, à égaliser les variances, et à transformer un effet multiplicatif en effet additif. Si vous ne voulez pas (ou ne pouvez pas) transformer vos données, vous pouvez utiliser l'équivalent non-paramétrique de l'ANOVA : le [[http://en.wikipedia.org/wiki/Kruskal–Wallis_one-way_analysis_of_variance|test de Kruskal-Wallis ]]. # Diagrammes de diagnostic opar <- par(mfrow=c(2,2)) plot(anov1) par(opar) # Test de la supposition de la normalité des résidus shapiro.test(resid(anov1)) # Test de la supposition de l'homogénéité de la variance bartlett.test(logMaxAbund ~ Diet, data=bird) Idéalement, le premier graphique devrait montrer une dispersion similaire pour chaque niveau de diète. Toutefois, les tests de Shapiro et de Bartlett ne sont pas significatifs. On peut supposer que les résidus sont distribués normalement et que les variances sont égales. ==== 3.5 Sortie du modèle ==== Lorsque le modèle d'ANOVA a été validé, on peut interpréter les résultats correctement. La sortie du modèle fournie par R dépend de la fonction qui a été utilisée pour effectuer l'ANOVA. Si la fonction ''aov()'' a été utilisée : aov1 <- aov(logMaxAbund ~ Diet, data=bird) Les résultats de l'ANOVA peuvent être visualisés avec la fonction ''summary()'' : summary(aov1) Si la fonction ''lm()'' a été utilisée : anov1 <- lm(logMaxAbund ~ Diet, data=bird) les résultats de l'ANOVA peuvent être visualisés avec la fonction ''anova()'' : anova(anov1) Dans les deux cas, la sortie dans R sera la même : Df Sum Sq Mean Sq F value Pr(>F) Diet 4 5.106 1.276 2.836 0.0341 * Residuals 49 22.052 0.450 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Cette sortie de R représente le tableau de l'ANOVA. On y retrouve les degrés de liberté, la somme des carrés, la moyenne de la somme des carrés, la statistique de F ainsi qu'une valeur de p. Dans l'exemple de la diète des oiseaux, la diète influence significativement l'abondance des oiseaux car la valeur de p est inférieure à 0.05. L'hypothèse nulle est rejetée, ce qui signifie qu'au moins une des diètes influence l'abondance différemment des autres diètes. ==== 3.6 Test complémentaire ==== Il est impossible d'identifier quel traitement diffère des autres avec une ANOVA. Elle ne permet que de déterminer s'il existe une différence entre niveaux. Pour identifier les niveaux qui diffèrent des autres, les tests post-hoc comparent les combinaisons de variables explicatives (i.e. les traitements) deux par deux. Il existe plusieurs tests post hoc (e.g. Fischer’s least significant difference, Duncan’s new multiple range test, Newman-Keuls method, Dunnett’s test, etc.), mais le test de Tukey est le plus couramment utilisé. Dans R, on utilise la fonction ''TukeyHSD()'' pour effectuer ce test : # À quel niveau se situe la différence de diète ? TukeyHSD(aov(anov1),ordered=T) # Cette commande est équivalente à la précédente : TukeyHSD(aov1,ordered=T) La sortie de R inclut un tableau qui liste toutes les combinaisons des niveaux de la variable explicative et qui identifie quel(s) traitement(s) diffère(ent) des autres :\\ Tukey multiple comparisons of means 95% family-wise confidence level factor levels have been ordered Fit: aov(formula = anov1) $Diet diff lwr upr p adj Vertebrate-InsectVert 0.3364295 -1.11457613 1.787435 0.9645742 Insect-InsectVert 0.6434334 -0.76550517 2.052372 0.6965047 Plant-InsectVert 0.8844338 -1.01537856 2.784246 0.6812494 PlantInsect-InsectVert 1.0657336 -0.35030287 2.481770 0.2235587 Insect-Vertebrate 0.3070039 -0.38670951 1.000717 0.7204249 Plant-Vertebrate 0.5480043 -0.90300137 1.999010 0.8211024 PlantInsect-Vertebrate 0.7293041 0.02128588 1.437322 0.0405485 Plant-Insect 0.2410004 -1.16793813 1.649939 0.9884504 PlantInsect-Insect 0.4223003 -0.19493574 1.039536 0.3117612 PlantInsect-Plant 0.1812999 -1.23473664 1.597336 0.9961844 Dans ce cas-ci, la seule différence significative d'abondance se retrouve entre les diètes "PlantInsect" et "Vertebrate". ==== 3.7 Visualisation des résultats==== Après avoir vérifié les conditions de base, interprété les résultats, et identifié les niveaux significatifs à l'aide de tests post-hoc ou de contrastes, les résultats d'une ANOVA peuvent être représentés graphiquement à l'aide de la fonction ''barplot()''. R produit donc un graphique de la variable réponse en fonction des niveaux de traitement, ou les erreurs types et le nom des niveaux du traitement (représentant le résultat d'un test post hoc) peuvent y être apposées. # Visualisation d'un modèle d'ANOVA à l'aide de la fonction barplot() sd <- tapply(bird$logMaxAbund,list(bird$Diet),sd) means <- tapply(bird$logMaxAbund,list(bird$Diet),mean) n <- length(bird$logMaxAbund) se <- 1.96*sd/sqrt(n) bp <- barplot(means, col=c("white","lightblue1","skyblue1","skyblue3","skyblue4"), ylab = expression("log"[10]*"(Maximum Abundance)"), xlab="Diet", ylim=c(0,1.8)) # Ajout des lignes verticales représentant les erreurs types segments(bp, means - se, bp, means + se, lwd=2) # et des lignes horizontales segments(bp - 0.1, means - se, bp + 0.1, means - se, lwd=2) segments(bp - 0.1, means + se, bp + 0.1, means + se, lwd=2) {{:lm_fig_12.png?650|}} ==== 3.8 Contrastes (section avancée/facultative) ==== * Les contrastes sont des comparaisons de moyennes basées sur des hypothèses //a priori//, * Ces groupes peuvent être composés d'un ou plusieurs niveaux d'un facteur, * On peut tester une hypothèse simple (e.g. μ₁ = μ₂) ou des hypothèses plus complexes (e.g. (μ₁ + μ₂)/3 == μ₃). Le nombre de comparaisons doit être plus petit ou égal au nombre de degrés de liberté de l'ANOVA. Ces comparaisons doivent être indépendantes l'une de l'autre. Il est possible d'afficher des résultats supplémentaires de l'ANOVA dans R, qu'on appelle contrastes. Ceci permet de visualiser les estimations de paramètres pour chaque niveau de la variable catégorique en comparaison avec un niveau de référence. On peut afficher ces résultats avec la fonction ''summmary.lm()'' lorsque l'ANOVA a été effectuée avec la fonction ''aov()'' et avec la fonction ''summary()'' lorsque l'ANOVA a été effectuée avec la fonction ''lm()''. Cette sortie montre les résultats de la régression linéaire pour chaque niveau de la variable catégorique : Call: lm(formula = logMaxAbund ~ Diet, data = bird) Residuals: Min 1Q Median 3Q Max -1.85286 -0.32972 -0.08808 0.47375 1.56075 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.1539 0.1500 7.692 5.66e-10 *** DietInsectVert -0.6434 0.4975 -1.293 0.2020 DietPlant 0.2410 0.4975 0.484 0.6303 DietPlantInsect 0.4223 0.2180 1.938 0.0585 . DietVertebrate -0.3070 0.2450 -1.253 0.2161 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6709 on 49 degrees of freedom Multiple R-squared: 0.188, Adjusted R-squared: 0.1217 F-statistic: 2.836 on 4 and 49 DF, p-value: 0.0341 La dernière ligne de cette sortie est identique à la sortie précédente de l'ANOVA. La statistique de F de l'ANOVA et sa valeur de p associée (2.836 et 0.0341 respectivement) sont les mêmes que celles présentées dans la table d'ANOVA, ce qui indique que la variabilité de l'abondance est mieux expliquée par la diète que par un modèle nul; la diète a donc un effet significatif sur l'abondance. L'ajustement du modèle (i.e. le R² ajusté) apparaît sur l'avant-dernière ligne de la sortie. Dans ce cas-ci, la diète explique 12.17% de la variabilité de l'abondance.\\ Les contrastes sont utilisés pour ajuster une variable réponse en fonction des différents niveaux d'une variable catégorique. Dans le cas de l'abondance en fonction de la diète, cinq régressions linéaires (correspondant aux cinq coefficients dans la sortie de R) sont calculées par la fonction ''lm()'', car la variable "diète" contient cinq niveaux. Par défaut, le niveau de référence correspond au premier niveau (en ordre alphabétique) de la variable catégorique. Ce niveau est indiqué par la ligne ''intercept'' dans la sortie de R, soit ''Insect'' dans ce cas-ci. Ici, le coefficient estimé du niveau de référence est comparé à 0 par un test de t, alors que les autres coefficients sont comparés à celui du niveau de référence. Dans ce cas-ci, seule la diète ''PlantInsect'' est différente de la diète ''Insect'' (valeur de p = 0.0585).\\ En d'autres mots, cette sortie de R permet de calculer la moyenne de la variable réponse pour chaque niveau de diète. Par exemple : LogMaxAbund = 1.1539 pour la diète ''Insect'',\\ LogMaxAbund = 1.1539 – 0.6434 pour la diète ''InsectVert'',\\ LogMaxAbund = 1.1539 + 0.2410 pour la diète ''Plant'',\\ etc. Ce type de contrastes compare chaque niveau de la variable explicative à un niveau de référence. Dans R, ceci correspond à la fonction ''contr.treatment()'' et représentent la méthode par défaut de la fonction ''lm()''. Le niveau de référence peut être changé en utilisant la fonction ''relevel()''. Par exemple, bird$Diet2 <- relevel(bird$Diet, ref="Plant") anov2 <- lm(logMaxAbund ~ Diet2, data=bird) summary(anov2) anova(anov2) compare chaque diète à la diète ''Plant'' maintenant définie comme le niveau de référence. La matrice de coefficients de contrastes peut être affichée par la commande suivante : contrasts(bird$Diet2) Insect InsectVert PlantInsect Vertebrate Plant 0 0 0 0 Insect 1 0 0 0 InsectVert 0 1 0 0 PlantInsect 0 0 1 0 Vertebrate 0 0 0 1 où chaque colonne représente une comparaison avec le niveau de référence ''Plant'' et chaque ligne représente un type de diète. Par exemple, la première comparaison est effectuée entre les diètes ''Plant'' et ''Insect''. La seconde comparaison est effectuée entre ''Plant'' et ''InsectVert'' etc. Il est possible de créer une matrice de coefficients de contrastes afin d'effectuer certaines comparaisons bien précises à l'aide de la fonction ''contrasts()''. Par exemple, contrasts(bird$Diet2) <- cbind(c(4,-1,-1,-1,-1), c(0,1,1,-1,-1), c(0,0,0,1,-1), c(0,1,-1,0,0)) crée la matrice de coefficients de contrastes suivante : [,1] [,2] [,3] [,4] Plant 4 0 0 0 Insect -1 1 0 1 InsectVert -1 1 0 -1 PlantInsect -1 -1 1 0 Vertebrate -1 -1 -1 0 qui compare :\\ * la diète ''Plant'' à toutes les autres diètes dans la première comparaison (colonne),\\ * les diètes ''InsectVert'' et ''Insect'' aux diètes ''PlantInsect'' et ''Vertebrate'' dans la deuxième comparaison,\\ * la diète ''PlantInsect'' à la diète ''Vertebrate'' dans la troisième comparaison,\\ * et la diète ''Insect'' à la diète ''InsectVert'' dans la quatrième comparaison. Pour chaque colonne, les diètes avec le même coefficient appartiennent au même groupe (e.g. pour la colonne 1, les quatre diètes avec un coefficient de -1 appartiennent au même groupe et sont comparés à la diète avec un coefficient différent; ici la diète "Plant" avec un coefficient de 4). Il est donc possible d'effectuer n'importe quelle comparaison possible à l'aide d'une matrice de coefficients de contrastes. Deux conditions doivent être respectées afin d'utiliser correctement ces matrices : - Pour chaque colonne, la somme des coefficients doit être égale à zéro et - la somme des produits de chaque paire de colonnes doit être égale à zéro. Ceci peut être vérifié à l'aide de la commande suivante : sum(contrasts(bird$Diet)[,1]) # première condition pour la colonne 1 sum(contrasts(bird$Diet)[,1]*contrasts(bird$Diet)[,2]) # deuxième condition pour les colonnes 1 et 2 Les contrastes utilisés fréquemment sont déjà programmés dans R : contrastes de Helmert, contrastes polynomiaux, etc. (cf. ''help(contrasts)''). ---- ===== 4. ANOVA à deux critères de classification ===== Jusqu'ici, les modèles ANOVA que nous avons explorés n'ont eu qu'une seule variable catégorique, mais on peut aussi créer des modèles ANOVA avec plusieurs variables explicatives catégoriques. Quand il y a 2 variables explicatives catégoriques, le modèle est un **ANOVA à deux critères de classification**. Un ANOVA à deux critères de classification teste plusieurs **hypothèses**: que la moyenne diffère entre les niveaux de la variable A, que la moyenne ne diffère pas entre les niveaux de la variable B; et qu'il n'y a pas d'interaction entre les variables A et B. Une **interaction significative** implique que la valeur moyenne de la variable réponse pour chaque niveau de la variable A change selon le niveau de la variable B. Par exemple, la relation entre la couleur d'un fruit et sa masse dépend de l'espèce de la plante: si oui, on peut dire qu'il y a une interaction entre la couleur et l'espèce. Cliquez ci-dessous pour plus des détails mathématiques. Le tableau de calcul de l'ANOVA vu à la section précédente doit être modifié afin d'inclure le deuxième facteur et l'interaction entre ces deux facteurs : ^ Source de\\ variation^ Degrés de\\ liberté (df)^ Sommes des carrés\\ des écarts à la moyenne^ Moyenne des carrés ^ Statistique de F^ ^ Total |abr-1 | {SS_Tot}=sum{i,j,k}{}({y_ijk}-overline{y})^2 | | | ^ Intra-\\ cases (erreur)|ab(r-1) | {SS_E}= sum{i,j,k}{}({y}_ijk-{overline{y}_ij})^2 |{MS_E}={SS_E}/{ab(r-1)}| | ^ Cases |ab-1 | {SS_Cells}= sum{i,j}{}({overline{y}_ij}-overline{y})^2 | | | ^ Facteur A |a-1 | {SS_FacteurA}= rb sum{i}{}({overline{y}_i.}-overline{y})^2|{MS_FacteurA}={SS_FacteurA}/{a-1}|{F_FacteurA}={MS_FacteurA}/{MS_E}| ^ Facteur B |b-1 | {SS_FacteurB}= ra sum{j}{}({overline{y}_.j}-overline{y})^2|{MS_FacteurB}={SS_FacteurB}/{b-1}|{F_FacteurB}={MS_FacteurB}/{MS_E}| ^ Interaction\\ entre A et B |(a-1)(b-1)| {SS_AB}= r sum{i,j,k}{}({overline{y}_{..k}}-{overline{y}_{.jk}}-{overline{y}_{i.k}})^2| {MS_AB}={SS_AB}/{(a-1)(b-1)}|{F_AB}={MS_AB}/{MS_E}| a: nombre de niveaux de la variable explicative A; b: nombre de niveaux de la variable explicative B; r: nombre de répétitions par traitement ==== 4.1 Effectuer une ANOVA à deux critères de classification ==== Dans R, une ANOVA à deux critères de classification est effectuée de la même manière qu'une ANOVA à un critère de classification avec la fonction ''lm()''. ---- **DÉFI 2** Examinez les effets des facteurs "Diet", "Aquatic" et de leur interaction sur l'abondance maximale d'oiseaux. Rappelez-vous que vous devez vérifier les suppositions statistiques de base avant d'interpréter les résultats d'une ANOVA, soit :\\ - Distribution normale des rsidus du modèle - Homoscédasticité des résidus de la variance Cette vérification peut être faite en utilisant les quatre graphiques de diagnostic expliqués dans la section précédente.\\ ++++ Défi 2: Solution | anov4 <- lm(logMaxAbund ~ Diet*Aquatic, data=bird) opar <- par(mfrow=c(2,2)) plot(anova4) par(opar) summary(anov4) anova(anov4) La fonction ''anova()'' permet de visualiser le tableau d'ANOVA du modèle :\\ Analysis of Variance Table Response: logMaxAbund Df Sum Sq Mean Sq F value Pr(>F) Diet 4 5.1059 1.27647 3.0378 0.02669 * Aquatic 1 0.3183 0.31834 0.7576 0.38870 Diet:Aquatic 3 2.8250 0.94167 2.2410 0.09644 . Residuals 45 18.9087 0.42019 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Dans ce cas-ci, le seul facteur significatif est la diète. La valeur de p de l'interaction n'est pas significative, ce qui signifie que l'effet de la diète est le même peu importe si l'oiseau est aquatique ou non. Le seuil de signification peut aussi être testé en comparant deux modèles nichés, i.e. en incluant un premier modèle avec une interaction et un deuxième modèle avec sans l'interaction. La fonction ''anova()'' est utilisée : anov5 <- lm(logMaxAbund ~ Diet + Aquatic, data=bird) anova(anov5, anov4) Voici la sortie dans R : Analysis of Variance Table Model 1: logMaxAbund ~ Diet + Aquatic Model 2: logMaxAbund ~ Diet * Aquatic Res.Df RSS Df Sum of Sq F Pr(>F) 1 48 21.734 2 45 1 8.909 3 2.825 2.241 0.09644 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Étant donné que la seule différence entre ces deux modèles est la présence de l'interaction, cette sortie de R présente le seuil de signification de cette interaction. Dans ce cas-ci, l'interaction n'est pas significative et peut donc être retirée du modèle.\\ Lorsque l'interaction est significative, rappelez-vous que chaque facteur ne peut pas être interprété séparément. Seule l'interaction peut l'être. __**Note :**__ Le tableau d'ANOVA indique que le nombre de degrés de liberté pour l'interaction entre la diète et le type d'oiseau (aquatique ou non) est de 3. Selon la notation mathématique du tableau de l'ANOVA à deux critères de classification (pour les plans équilibrés), a = 5 et b = 2 et le nombre de degrés de liberté est de (a-1)(b-1) = 4*1 = 4. La sortie de R indique cependant que le nombre de degrés de liberté est de 3. Cette interaction est extrêmement non équilibrée () : les oiseaux aquatiques ne se nourrissent pas de plante, donc ce niveau n'est pas considéré (prenez note du NA dans la sortie de summary(anov4)). Consultez la section avancée sur les [[http://qcbs.ca/wiki/r_atelier4#anova_non_equilibree_section_avancee_et_optionnelle|ANOVA non équilibrées]] plus bas pour plus de détails. ++++ ---- ==== 4.2 Diagramme d'interaction ==== Les interactions peuvent être visualisées à l'aide de la fonction ''interaction.plot()'' : interaction.plot(bird$Diet, bird$Aquatic, bird$logMaxAbund, col="black", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab="Diet") {{:interaction_plot.png?650|}} Que signifie le trou sur la ligne des oiseaux aquatiques? table(bird$Diet, bird$Aquatic) 0 1 Insect 14 6 InsectVert 1 1 Plant 2 0 PlantInsect 17 1 Vertebrate 5 7 Le plan est **non-équilibré**: il y a un nombre inégal d'observations entre les diètes pour les oiseaux aquatiques (représentés par le 1) et les oiseaux terrestres (représentés par le 0). Consultez la section avancée ci-dessous pour plus de détails sur les ANOVA à plan non-équilibré. ---- **DÉFI 3** Tester le seuil de signification du facteur "Aquatic" en comparant deux modèles nichés (i.e. avec et sans ce facteur). ++++ Défi 3: Solution | anova(anov3,anov5) # Rappel que anov3 est le modèle avec le facteur "Diet" seulement. ++++ ---- ===== 5. ANOVA non-équilibrée (section avancée/facultative) ===== Les ANOVA à un et à deux critères de classification nous permettent de déterminer les effets de variables catégoriques sur une variable réponse continue lorsque nous avons un plan expérimental **équilibré** (i.e. lorsque le nombre de répétitions est égal pour chaque niveau de facteur). Cependant, le plan peut devenir **non-équilibré** en raison d'une perte d'unités expérimentales au cours d'une expérience ou de restrictions techniques dues au plan expérimental. Dans de telles situations, les résultats d'ANOVA peuvent être mal interprétés en raison de calculs erronés de la somme des carrés. Pour des plans expérimentaux non-équilibrés, l'ANOVA doit être modifiée pour prendre en compte les données manquantes de la variable réponse. Le modèle mathématique, les hypothèses statistiques, et les conditions de base d'une ANOVA à plan non-équilibré demeurent les mêmes que deux de l'ANOVA à plan équilibré. Le calcul de la somme des carrés, par contre, est différent. Pour un plan expérimental non-équilibré, les hypothèses statistiques sont les suivantes : H₀: µ₁ = µ₂ =... = µ_i =... = µ_n\\ H₁: il y a au moins une moyenne µi qui diffère des autres. On utilise ce modèle mathématique : {y_{ijk}} = µ + {A_i} + {B_j} + {A_i}{B_j} + {ε_{ijk}} Rappelez-vous du calcul de la somme des carrés dans le cas d'une ANOVA à plan équilibré : {SS_FactorA} = rb sum{i}{}{({overline{y}_{i.}}-overline{y})}^2 = SS(A) \\ {SS_FactorB} = ra sum{j}{}{({overline{y}_{.j}}-overline{y})}^2 = SS(B|A) = SS(A,B)-SS(B) \\ {SS_AB} = r sum{i,j,k}{}({overline{y}_{..k}}-{overline{y}_{.jk}}-{overline{y}_{i.k}})^2 = SS(A,B,AB)-SS(A,B) Ceci correspond à une somme des carrés séquentielle (aussi appelée de type I), car l'effet d'un facteur B est calculé après avoir retiré l'effet d'un facteur A. L'interaction est calculée après avoir retiré les effets principaux de ces deux facteurs. Ces calculs dépendent de la taille de l'échantillon, car l'effet de chaque facteur est calculé après avoir retiré l'effet du facteur précédent. Dans le cas d'un plan expérimental non-équilibré, les résultats d'ANOVA dépendent de l'ordre dans lequel chaque variable apparaît dans le modèle. Voyez comment les résultats diffèrent en comparant les deux modèles suivants : unb_anov1 <- lm(logMaxAbund ~ Aquatic + Diet, data=bird) unb_anov2 <- lm(logMaxAbund ~ Diet + Aquatic, data=bird) anova(unb_anov1) anova(unb_anov2) Bien que les variables explicatives soient les mêmes pour chaque modèle, les tableaux d'ANOVA montrent des résultats différents à cause d'un plan non-équilibré (i.e. un nombre différent d'observations pour les oiseaux aquatiques et terrestres). Pour les plans non-équilibrés, une somme des carrés marginale (aussi appelée de type III) permet de calculer un effet principal après avoir retiré les effets des autres facteurs. Ceci rend le calcul indépendant de la taille des échantillons : {SS_FactorA}={SS(A|B,AB)}=SS(A,B,AB)-SS(B,AB)\\ {SS_FactorB}={SS(B|A,AB)}=SS(A,B,AB)-SS(A,AB)\\ {SS_AB}={SS(AB|B,A)}=SS(A,B,AB)-SS(B,AB) \\ \\ Dans R, une ANOVA avec somme des carrés de type III peut être effectuée avec la fonction ''Anova()'' du paquet "car" et en spécifiant l'argument ''type="III"'' : Anova(unb_anov1,type="III") En comparant les tableaux d'ANOVA de modèles avec un ordre différent dans les variables explicatives, on peut voir que les résultats sont les mêmes. L'utilisation de la somme des carrés de type III produit une ANOVA qui est indépendante de la tailles des échantillons. Après avoir vérifié les suppositions du modèle, les résultats peuvent finalement être interprétés correctement. ---- ===== 6. ANCOVA ===== L'analyse de covariance (ANCOVA) est un modèle linéaire qui teste l'influence d'une (ou plusieurs) variable explicative catégorique sur une variable réponse continue. Chaque niveau de la variable catégorique est décrit par une pente (ou coefficient de variation) et une ordonnée à l'origine. En plus de tester si la variable réponse diffère pour au moins un niveau de la variable catégorique, l'ANCOVA teste aussi si la variable réponse est influencée par sa relation avec la variable continue (nommée la **covariable** dans une ANCOVA), et par une différence dans l'influence de la variable continue sur la réponse (i.e. l'interaction) entre les niveaux de groupe. Les hypothèses d'un ANCOVA sont alors: qu'il n'y a pas de différence de moyenne entre les niveaux de la variable catégorique; qu'il n'y a pas de correlation entre la variable réponse et la variable explicative catégorique; et qu'il n'y a pas d'interaction entre les variables catégoriques et continues. ==== 6.1 Conditions de base ==== Tout comme le test de t et l'ANOVA, l'ANCOVA doit respecter certaines conditions statistiques qu'on peut vérifier à l'aide de diagrammes de diagnostic: - Les covariables ont toutes la **même étendue de valeurs** - Les variables sont **//fixes//** - Les variables catégoriques et continues sont **indépendantes** **__Note :__** Un variable //fixe// est une variable d'intérêt pour une étude (e.g. la masse des oiseaux). En comparaison, une variable aléatoire représente surtout une source de bruit qu'on veut contrôler (i.e. le site où les oiseaux ont été échantillonnés). Si votre modèle comporte des effets aléatoires, consultez l'atelier sur les [[http://qcbs.ca/wiki/r_atelier6|modèles linéaires mixtes]]! ==== 6.2 Types d'ANCOVA ==== Il est possible d'avoir plusieurs facteurs (i.e. variables explicatives catégoriques) et covariables (i.e. variables explicatives continues) au sein d'une même ANCOVA. Par contre, l'interprétation des résultats devient de plus en plus complexe à mesure que le nombre de covariables et de facteurs augmente. Les ANCOVA les plus courantes comportent : - une covariable et un facteur - une covariable et deux facteurs - deux covariables et un facteur Les buts possibles de l'ANCOVA sont de déterminer les effets : - des facteurs et des covariables sur la variable réponse - des facteurs sur la variable réponse après avoir retiré l'effet des covariables - des facteurs sur la relation existant entre les covariables et la variable réponse Ces buts ne sont atteints que s'il n'y a pas d'interaction significative entre le(s) facteur(s) et la(les) covariable(s)! Des exemples d'interaction significative entre un facteur et une covariable (pour une ANCOVA avec un facteur et une covariable) sont illustrés ci-dessous dans les deux derniers graphiques: {{:ancova_schematic.png?550|}} La même logique s'applique aux ANCOVAs à plusieurs facteurs et/ou covariables. \\ ==== 6.3 Effectuer une ANCOVA ==== Effectuer une ANCOVA dans R ressemble à une ANOVA à deux critères de classification : on utilise la fonction ''lm()''. Toutefois, au lieu d'avoir deux variables catégoriques (e.g. "Diet" et "Aquatic"), on utilise maintenant une variable catégorique et une variable continue. Par exemple, en utilisant le jeu de données CO2 (déjà inclus dans R) où la variable réponse est //uptake//, on peut effectuer une ANCOVA avec la variable continue //conc// et le facteur //Treatment// : ancova.example <- lm(uptake ~ conc*Treatment, data=CO2) anova(ancova.example) Si l'analyse indique que seule la covariable est significative, on retire le facteur du modèle; on revient à une **ANOVA** à un critère de classification.\\ Si l'analyse indique que seul le facteur est significatif, on retire la covariable du modèle; on revient à une **régression linéaire simple**.\\ Si l'analyse indique que l'interaction est significative, il faut trouver quels niveaux ont une pente différente.\\ Dans l'exemple du jeu de données CO2, la covariable et le facteur sont significatifs, mais l'interaction n'est pas significative. Si on remplace le facteur //Treatment// par le facteur //Type//, l'interaction devient significative. Si vous voulez comparer les moyennes de la variable réponse entre les facteurs, vous pouvez utiliser les moyennes ajustées qui sont calculées comme dans l'équation de l'ANCOVA et en tenant compte de l'effet de la covariable : install.packages("effects") library(effects) adj.means <- effect('Treatment', ancova.example) plot(adj.means) adj.means <- effect('conc*Treatment', ancova.example) plot(adj.means) ---- **DÉFI 4** Effectuez une ANCOVA qui teste l'effet du facteur //Diet//, de la covariable //Mass//, et de leur interaction sur la variable réponse //MaxAbund//.\\ ++++ Défi 4: Solution | # Si vous avez complété la section avancée sur les contrastes, vous devrez réinitialiser les contrastes # à l'aide de la fonction ''options()'' # Si vous n'avez pas complété la section avancée sur les contrastes, ignorez la première ligne du script. options(contrasts=c("contr.treatment", "contr.poly")) ancov1 <- lm(logMaxAbund ~ logMass*Diet, data=bird) summary(ancov1) anova(ancov1) R fournit la sortie suivante pour cette ANCOVA : Analysis of Variance Table Response: logMaxAbund Df Sum Sq Mean Sq F value Pr(>F) logMass 1 1.9736 1.97357 4.6054 0.03743 * Diet 4 3.3477 0.83691 1.9530 0.11850 logMass:Diet 4 2.9811 0.74527 1.7391 0.15849 Residuals 44 18.8556 0.42854 Dans ce cas-ci, l'interaction n'est pas significative, ce qui signifie que l'effet de la masse sur l'abondance maximale est le même peu importe la diète. L'interaction est retirée du modèle et l'ANCOVA devient : ancov2 <- lm(logMaxAbund ~ logMass + Diet, data=bird) R nous indique que la diète n'est pas significative non plus, donc ce terme est retiré du modèle. Notre modèle final devient donc une régression linéaire simple : lm2 <- lm(logMaxAbund ~ logMass, data=bird) Les résultats de l'analyse peuvent être représentés graphiquement. On trace la variable réponse en fonction de la variable explicative continue avec des points et des lignes de différentes couleurs pour les différents niveaux de la variable catégorique. Nous pouvons aussi tracez un diagramme représentant les pentes et ordonnées à l'origine d'une ANCOVA (le modèle ancov1 plus haut) à l'aide des fonctions ''abline()'' et ''coef()''. coef(ancov1) plot(logMaxAbund~logMass, data=bird, col=Diet, pch=19, ylab=expression("log"[10]*"(Maximum Abundance)"), xlab=expression("log"[10]*"(Mass)")) abline(a=coef(ancov1)[1],b=coef(ancov1)[2], col="deepskyblue1") abline(a=sum(coef(ancov1)[1]+coef(ancov1)[3]),b=sum(coef(ancov1)[2]+coef(ancov1)[7]),col="green2", lwd=2) abline(a=sum(coef(ancov1)[1]+coef(ancov1)[4]),b=sum(coef(ancov1)[2]+coef(ancov1)[8]),col="orange1", lwd=2) abline(a=sum(coef(ancov1)[1]+coef(ancov1)[5]),b=sum(coef(ancov1)[2]+coef(ancov1)[9]),col="lightsteelblue1", lwd=2) abline(a=sum(coef(ancov1)[1]+coef(ancov1)[6]),b=sum(coef(ancov1)[2]+coef(ancov1)[10]),col="darkcyan", lwd=2) {{:ancova_plot.png?350|}} ++++ ---- ===== 7. Régression multiple ===== Une régression multiple teste les effets de plusieurs variables explicatives continues sur une variable réponse continue. ==== 7.1 Conditions de base ==== En plus des conditions de base habituelles des modèles linéaires, il est important de tester pour **l'orthogonalité** parce que ceci influence l'interprétation du modèle. Les variables sont orthogonales quand les variables explicatives sont colinéaires. Si une variable explicative est corrélée avec une autre variable, elles expliquent probablement la même variabilité dans la variable réponse, et l'effet d'une variable sera caché par l'autre.S'il existe une relation entre deux variables explicatives, elles sont donc **colinéaires**. La colinéarité doit être évitée, car il ne sera pas possible de distinguer les effets propres à chaque variable. Voici quelques solutions possibles : - Gardez seulement **une** des variables colinéaires, - Essayez une analyse multidimensionelle //[[http://qcbs.ca/wiki/r_atelier9|(voir l'atelier 9)]]//, - Essayez une analyse pseudo-orthogonale. La colinéarité entre variables explicatives peut être vérifiée à l'aide du **facteur d'inflation de la variance** (VIF ou //"variance inflation factor"//) en utilisant la fonction ''vif'' du paquet ''HH''. vif(clDD ~ clFD + clTmi + clTma + clP + grass, data=Dickcissel) qui produit la sortie suivante: clFD clTmi clTma clP grass 13.605855 9.566169 4.811837 3.196599 1.165775 Les variables colinéaires ont un VIF plus grand que 5. La sortie R montre alors que //clDD//, //clFD//, et //clTmi// sont fortement colinéaires. Seulement une de ces variables explicatives peut être retenue dans le modèle de régression final. ==== 7.2 Jeu de données Dickcissel ==== Le jeu de données Dickcissel (du nom d'un petit oiseau granivore de la famille des Cardinalidae) explore les effets de plusieurs variables environnementales qui pourraient expliquer l'abondance et la présence d'une espèce d'oiseau des prairies nord-américaines avec des pics d'abondance au Kansas, É-U. Le jeu de données contient 15 variables : ^ Nom de la variable^ Description^ Type^ ^ abund | Le nombre d'individus\\ observé sur chaque route | Continu/ numérique| ^ Present | Présence/ absence\\ de l'espèce| Binaire ("Présent"/ "Absent")| ^ broadleaf, conif, crop, grass,\\ shrub, urban, wetland | Variables du paysage à moins de 20 km de rayon\\ du centre de la route|Continu/ numérique| ^ NDVI | Indice de végétation (une mesure de la productivité)| Nombre entier| ^ clDD, clFD, clTma, clTmi, clP | Données climatiques (DD = degrés jours,\\ FD = jours de gel, Tma = température max,\\ Tmi = température min,\\ P = précipitation)| Continu/ numérique| Dans R, une régression multiple est effectuée à l'aide de la fonction ''lm()'' et les résultats peuvent être visualisés avec la fonction ''summary()''. En utilisant le jeu de données Dickcissel, on peut tester les effets du //climat//, de la //productivité// et du //paysage// sur l'//abondance// du Dickcissel à l'aide d'un modèle de régression multiple. ---- **DÉFI 5** Est-il nécessaire de transformer la variable réponse //abund// ? ++++ Défi 5: Solution| hist(Dickcissel$abund, main="", xlab="Dickcissel abundance") shapiro.test(Dickcissel$abund) skewness(Dickcissel$abund) summary(Dickcissel$abund) Il y a plusieurs zéros dans la distribution de la variable //abund//. On peut ajouter une constante avant d'effectuer une transformation logarithmique, étant donnée la forte asymétrie de la distribution : hist(log10(Dickcissel$abund+0.1), =main="", xlab=expression("log"[10]*"(Dickcissel Abundance + 0.1)")) shapiro.test(log10(Dickcissel$abund+0.1)) skewness(log10(Dickcissel$abund+0.1)) La variable n'est toujours pas distribuée normalement après la transformation. ++++ ---- Dans le défi 5, vous avez probablement remarqué que la variable réponse //abund// n'a pas pu être normalisée. Il faudrait alors peut-être laisser tomber la supposition de normalité et passer à un [[http://qcbs.ca/wiki/r_atelier7| modèle linéaire généralisé]], mais ceci ira à plus tard ! Pour l'instant, nous allons utiliser la variable //abund// non-transformée et comparer l'importance relative de trois variables explicatives (//climat//, //productivité// et //paysage//) sur l'//abondance//. lm.mult <- lm(abund ~ clTma + NDVI + grass, data=Dickcissel) summary(lm.mult) La sortie de R indique quelles variables explicatives sont significatives : lm(formula = abund ~ clTma + NDVI + grass, data = Dickcissel) Residuals: Min 1Q Median 3Q Max -35.327 -11.029 -4.337 2.150 180.725 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -83.60813 11.57745 -7.222 1.46e-12 *** clTma 3.27299 0.40677 8.046 4.14e-15 *** NDVI 0.13716 0.05486 2.500 0.0127 * grass 10.41435 4.68962 2.221 0.0267 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 22.58 on 642 degrees of freedom Multiple R-squared: 0.117, Adjusted R-squared: 0.1128 F-statistic: 28.35 on 3 and 642 DF, p-value: < 2.2e-16 Dans ce cas-ci, les trois variables explicatives ont une influence significative sur l'abondance de Dickcissel, la plus significative étant le climat (p = 4.14e-15). Ces trois variables expliquent 11.28% de la variabilité de l'abondance de Dickcissel (R carré ajusté = 0.1128). Le modèle global est également significatif et explique la variabilité de l'abondance de Dickcissel mieux qu'un modèle nul (p < 2.2e-16). Un graphique de la variable réponse en fonction de chaque variable explicative peut être utilisé pour représenter les résultats du modèle : plot(abund ~ clTma, data=Dickcissel, pch=19, col="orange") plot(abund ~ NDVI, data=Dickcissel, pch=19, col="skyblue") plot(abund ~ grass, data=Dickcissel, pch=19, col="green") {{:lm_fig 13.png?850|}} ==== 7.3 Régression polynomiale (section avancée/facultative) ==== Les relations entre variables réponses et variables explicatives ne sont pas toujours linéaires. Dans ce cas, une régression linéaire, qui correspond à une droite reliant deux variables, est incapable de représenter la relation entre ces deux variables. Une régression polynomiale permet de tracer une courbe polynomiale entre la variable réponse et la variable explicative, et permet de représenter une relation non linéaire basée sur le modèle mathématique suivant : {y_i} = {β_0} + {β_1}{x_i} + {β_2}{{x_i}^2} + {β_3}{{x_i}^3} + {ε_i} pour un polynôme de degré 3\\ {y_i} = {β_0} + {β_1}{x_i} + {β_2}{{x_i}^2} + {ε_i} pour un polynôme de degré 2 où β₀ est l'ordonnée à l'origine de la droite de régression,\\ β₁ est l'effet de la variable x,\\ β₂ est l'effet de la variable x au carré (x²),\\ β₃ est l'effet de la variable x au cube (x³),\\ ε_i sont les résidus du modèle (i.e. la variation inexpliquée). Le degré du polynôme est l'exposant le plus élevé de l'équation. En connaissant le degré du polynôme, on peut le qualifier : {{:polynomial_degree.png?250|}} Dans R, ces modèles sont effectués avec la fonction ''lm()'' et peuvent être comparés entre eux avec la fonction ''anova()'' : lm.linear <- lm(abund ~ clDD, data=Dickcissel) lm.quad <- lm(abund ~ clDD + I(clDD^2), data=Dickcissel) lm.cubic <- lm(abund ~ clDD + I(clDD^2) + I(clDD^3), data=Dickcissel) ---- **DÉFI 7** Comparez les différents modèles polynomiaux de l'exemple précédent et déterminez quel modèle est le plus approprié. Extrayez le R carré ajusté, les coefficients de régression et les valeurs de p de ce modèle. ++++ Défi 7: Solution| anova(lm.linear,lm.quad,lm.cubic) # fait la liste des modèles par ordre croissant de complexité. # On accepte le modèle se trouvant sur la ligne la plus basse avec une valeur de p significative. # i.e. le modèle lm.quad # Examinons le résumé du modèle summary(lm.quad) # Les coefficients de régression summary(lm.quad)$coefficients[,1] # Estimation des valeurs de p summary(lm.quad)$coefficients[,4] # Le R carré ajusté summary(lm.quad)$adj.r.squared ++++ ---- La comparaison de modèles du défi 7 a montré que la régression quadratique (i.e. polynomiale de degré 2) était le meilleur modèle. Le polynôme de degré trois peut être retiré du modèle final : Analysis of Variance Table Model 1: abund ~ clDD Model 2: abund ~ clDD + I(clDD^2) Model 3: abund ~ clDD + I(clDD^2) + I(clDD^3) Model Res.Df RSS Df Sum of Sq F Pr(>F) 1 644 365039 2 643 355871 1 9168.3 16.5457 5.34e-05 *** 3 642 355743 1 127.7 0.2304 0.6314 La sortie de R pour le modèle final est : Call: lm(formula = abund ~ clDD + I(clDD^2), data = Dickcissel) Residuals: Min 1Q Median 3Q Max -14.057 -12.253 -8.674 1.495 190.129 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.968e+01 5.954e+00 -3.306 0.001 ** clDD 1.297e-02 2.788e-03 4.651 4.00e-06 *** I(clDD^2) -1.246e-06 3.061e-07 -4.070 5.28e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 23.53 on 643 degrees of freedom Multiple R-squared: 0.04018, Adjusted R-squared: 0.0372 F-statistic: 13.46 on 2 and 643 DF, p-value: 1.876e-06 Dans cet exemple, le terme linéaire influence plus la variable réponse que le terme quadratique, car leurs valeurs de p sont 4.00e-06 et 5.28e-05 respectivement. Ces deux termes expliquent 3.72% de la variabilité de l'abondance (R carré ajusté), ce qui est très peu. ==== 7.4 Régression pas à pas ==== Afin d'obtenir un modèle de régression multiple "optimal", on peut commencer avec un modèle qui inclut toutes les variables explicatives et retirer les variables non significatives, procédant donc à une **sélection pas à pas**. Les variables non significatives sont retirées une à la fois et l'ajustement de chaque modèle successif est évalué à l'aide de l'AIC [[http://fr.wikipedia.org/wiki/Crit%C3%A8re_d%27information_d%27Akaike|(Critère d'information d'Akaike)]], jusqu'à ce que toutes les variables explicatives soient significatives. Prenez note qu'une valeur plus basse d'AIC indique un meilleur ajustement (i.e. le meilleur modèle est celui avec la valeur d'AIC la plus basse). Dans R, la sélection pas à pas est effectuée à l'aide de la fonction ''step()'' : lm.full <- lm(abund ~ . - Present, data=Dickcissel) lm.step <- step(lm.full) summary(lm.step) Les résultats indiquent que seulement 6 variables sont significatives parmi les 13 variables de départ : Call: lm(formula = abund ~ clDD + clFD + clTmi + clTma + clP + grass, data = Dickcissel) Residuals: Min 1Q Median 3Q Max -30.913 -9.640 -3.070 4.217 172.133 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -4.457e+02 3.464e+01 -12.868 < 2e-16 *** clDD 5.534e-02 8.795e-03 6.292 5.83e-10 *** clFD 1.090e+00 1.690e-01 6.452 2.19e-10 *** clTmi -6.717e+00 7.192e-01 -9.339 < 2e-16 *** clTma 3.172e+00 1.288e+00 2.463 0.014030 * clP 1.562e-01 4.707e-02 3.318 0.000959 *** grass 1.066e+01 4.280e+00 2.490 0.013027 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 19.85 on 639 degrees of freedom Multiple R-squared: 0.3207, Adjusted R-squared: 0.3144 F-statistic: 50.29 on 6 and 639 DF, p-value: < 2.2e-16 Le modèle explique maintenant 31.44% de la variabilité de l'abondance. La variable explicative la plus significative est clTmi.\\ Cependant, certaines des variables explicatives sont corrélées entre elles et devraient être retirées du modèle final afin de ne pas inclure de variables qui n'apportent pas d'information nouvelle. ===== 8. Partition de la variation (section avancée/facultative) ===== Afin d'évaluer la contribution relative de deux ou plusieurs variables explicatives à la variabilité d'une variable réponse, on peut utiliser la fonction ''varpart()'' de la librairie "vegan". Cette fonction permet de subdiviser la variation expliquée de la réponse variable entre différents groupes de variables explicatives. Par exemple, dans le jeu de données Dickcissel, on peut évaluer les contributions relatives des données climatiques et du paysage de la manière suivante : part.lm = varpart(Dickcissel$abund, Dickcissel[,c("clDD", "clFD", "clTmi", "clTma", "clP")], Dickcissel[,c("broadleaf", "conif", "grass", "crop", "urban", "wetland")]) part.lm La sortie de R permet de visualiser la partition de la variation : Partition of variation in RDA Call: varpart(Y = Dickcissel$abund, X = Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")], Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")]) Explanatory tables: X1: Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")] X2: Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")] No. of explanatory tables: 2 Total variation (SS): 370770 Variance: 574.84 No. of observations: 646 Partition table: Df R.squared Adj.R.squared Testable [a+b] = X1 5 0.31414 0.30878 TRUE [b+c] = X2 6 0.03654 0.02749 TRUE [a+b+c] = X1+X2 11 0.32378 0.31205 TRUE Individual fractions [a] = X1|X2 5 0.28456 TRUE [b] 0 0.02423 FALSE [c] = X2|X1 6 0.00327 TRUE [d] = Residuals 0.68795 FALSE --- Utilisez la fonction ''rda'' pour tester si les différentes fractions sont significatives. Cette sortie R montre que les deux groupes de variables explicatives expliquent 31.205% ([a+b+c] = X1+X2) de la variation de l'abondance de Dickcissel alors que les variables du climat expliquent à elles seules 28.46% de la variation ([a] = X1|X2) et les variables du paysage ne contribuent qu'à expliquer 0.33% de la variation de l'abondance de Dickcissel ([c] = X2|X1). L'interaction entre les deux groupes de variables expliquent 2.42% ([b]) de la variation. Pour tester si chaque fraction est significative, il est possible d'utiliser la RDA partielle et un test par permutation avec les fonctions ''rda()'' et ''anova()'' : # Variables climatiques out.1 = rda(Dickcissel$abund, Dickcissel[,c("clDD", "clFD", "clTmi", "clTma", "clP")], Dickcissel[,c("broadleaf", "conif", "grass", "crop", "urban", "wetland")]) anova(out.1, step=1000, perm.max=1000) # Variables du paysage out.2 = rda(Dickcissel$abund, Dickcissel[,c("broadleaf", "conif", "grass", "crop", "urban", "wetland")], Dickcissel[,c("clDD", "clFD", "clTmi", "clTma", "clP")]) anova(out.2, step=1000, perm.max=1000) la sortie R : Variables climatiques Permutation test for rda under reduced model Model: rda(X = Dickcissel$abund, Y = Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")], Z = Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")]) Df Var F N.Perm Pr(>F) Model 5 165.12 53.862 999 0.001 *** Residual 634 388.72 Variables du paysage Permutation test for rda under reduced model Model: rda(X = Dickcissel$abund, Y=Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")], Z = Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")]) Df Var F N.Perm Pr(>F) Model 6 5.54 1.5063 999 0.152 Residual 634 388.72 Dans ce cas, la fraction de la variation expliquée par les variables du climat est significative (p-value=0.001) alors que la fraction expliquée par les variables du paysage ne l'est pas (p-value=0.152). Les résultats du partitionnement de la variation sont généralement représentés graphiquement par un diagramme de Venn dans lequel chaque groupe de variables explicatives est représenté par un cercle. La fraction de la variation expliquée est indiquée à l'intérieur des cercles. showvarparts(2) plot(part.lm,digits=2, bg=rgb(48,225,210,80,maxColorValue=225), col="turquoise4") {{::varpart.png?300|}} ---- ===== Allez plus loin ! ===== **Super ! Vous êtes maintenant prêts à effectuer des régressions, des ANOVA et des ANCOVA sur vos propres données. Cependant, rappelez-vous de toujours spécifier vos modèles correctement et de vérifier leurs conditions de base avant d'interpréter les résultats en fonction des caractéristiques écologiques de vos données. ** Quelques livres pertinents à propos des régressions linéaires et de l'ANOVA * Myers RH - Classical and Modern Regression with Application * Gotelli NJ - A Primer of Ecological Statistics