This is an old revision of the document!
Ateliers R du CSBQ
Cette série de 10 ateliers guide les participants à travers les étapes requises afin de maîtriser le logiciel R pour une grande variété d’analyses statistiques pertinentes en recherche en biologie et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R.
Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil
Atelier 4 : Modèles linéaires
Développé par : Catherine Baltazar, Bérenger Bourgeois, Zofia Taranu
Résumé : Durant 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 présentation Prezi associée : Prezi
Téléchargez les scripts R et les données pour cet atelier :
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
- Partition de la variance
(section avancée et optionnelle)
Aperçu
Les scientifiques sont souvent intéressés à déterminer les relations entre des variables. Selon la nature et le nombre de variables considérées, différents outils statistiques peuvent être utilisés pour évaluer ces relations. La table suivante dresse une liste de cinq types d'analyses statistiques qui seront couverts durant 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. Régression linéaire simple
Le but de cette analyse est de trouver une relation entre une variable réponse y et une variable explicative x en faisant passer une droite entre les données. Le modèle mathématique correspondant à une régression linéaire est représenté par l'équation suivante :
où
est l'ordonnée à l'origine de la droite de régression,
est la pente de la droite de régression,
est la variable explicative continue pour la iième observation,
sont les résidus du modèle (i.e. la variance inexpliquée).
L'objectif est de trouver la meilleure estimation de ces deux paramètres de régression (i.e. la pente et l'ordonnée à l'origine) et d'évaluer l'ajustement (i.e. goodness of fit) du modèle de régression. Bien que plusieurs méthodes aient été développées afin de calculer les coefficients de la pente et de l'ordonnée à l'origine d'un modèle de régression, la méthode des moindres carrés est la méthode la plus utilisée et correspond à la méthode par défaut de 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. La pente (β1) et l'ordonnée à l'origine (β0) peuvent être calculées de la façon suivante :
Afin d'être valide, une régression linéaire doit respecter quatre suppositions de base sans lesquelles le modèle ne peut pas être interprété correctement.
1.1 Suppositions de base
- Homoscédasticité
Les variables explicatives doivent avoir une variance homogène (également appelée homoscédasticité), i.e. la dispersion des données doit être uniforme pour chaque valeur de xi. Cette supposition peut être vérifiée indirectement en représentant graphiquement la dispersion des résidus en fonction des valeurs prédites de la variable réponse.
Nous verrons plus loin qu'en cas d'hétéroscédasticité, une transformation des données ou un modèle linéaire généralisé avec une distribution différente (Poisson, binomiale négative, etc.) peuvent être appliqués afin de mieux traduire la relation entre les variables. - Indépendance
Une régression linéaire peut seulement être appliquée à des données indépendantes. Ceci signifie qu'une valeur de yi pour une valeur de xi donnée ne doit pas être influencée par d'autres valeurs de xi. Le non-respect de cette supposition peut se produire lorsqu'il existe une forme de dépendance au sein des données telle qu'une corrélation spatiale ou temporelle. - Influence
Si des observations dans le jeu de données sont très différentes des autres observations, il peut y avoir des problèmes lors de la calibration du modèle, car ces observations vont fortement influencer la valeur de la pente et de l'ordonnée à l'origine. - Distribution normale
Une régression linéaire devrait seulement être appliquée à des données suivant une distribution normale (variables réponse et explicative).
1.2 Effectuer un modèle linéaire
Nous allons effectuer une première régression linéaire en analysant la relation entre l'abondance maximale et la masse du jeu de données “bird”.
Dans R, une régression linéaire est codé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. |
- | Charger et explorer les données
# Chargez les librairies 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 | L'abondance la plus élevée observée à n'importe quel site en Amérique du Nord | Continue/numérique |
AvgAbund | L'abondance moyenne sur tous les sites en Amérique du Nord | Continue/numérique |
Mass | La 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) |
Nous sommes maintenant prêts à exécuter le modèle linéaire :
- | Régression de l'abondance maximale en fonction de la masse
lm1 <- lm(bird$MaxAbund ~ bird$Mass) # où Y ~ X signifie Y "en fonction de" X>
1.3 Vérification des suppositions
- | Graphiques de diagnostic
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
Homoscédasticité
Graphique des résidus en fonction des valeurs prédites - Le premier graphique de diagnostic (créé avec la fonction plot(lm1)
) représente la dispersion des résidus en fonction des valeurs prédites par le modèle de régression linéaire. Ceci permet de vérifier si la condition d'homoscédasticité est respectée : ce graphique devrait montrer une dispersion similaire le long des valeurs prédites (axe des x). Si la relation entre la variable réponse et la variable explicative n'est pas linéaire, ce graphique va nous l'indiquer.
Graphique “Scale-location” - Le troisième graphique permet de vérifier si la dispersion des résidus augmente pour une valeur prédite donnée (i.e. ça identife si la dispersion des résidus est causée par la variable explicative). Si la dispersion augmente, la supposition d'homoscédasticité n'est pas respectée.
Indépendance et distribution normale
Diagramme quantile-quantile - L'indépendance peut être évaluée à l'aide d'un diagramme quantile-quantile. Ceci permet de vérifier la distribution des résidus du modèle et de vérifier la normalité de la variable réponse. 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 sont situés près d'une ligne 1:1, les résidus peuvent être considérés comme normalement distribués.
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
Résidus vs diagramme d'influence - L'influence de certaines données peut être visualisée sur le quatrième graphique (i.e. résidus en fonction de l'influence) qui identifie les numéros d'observations avec une haute influence. Si (et seulement si!) ces observations correspondent à des erreurs de mesure ou à des exceptions, elles peuvent être retirées du jeu de données.
1.4 La normalisation des données
Dans l'exemple précédent, la variable réponse MaxAbund et la variable explicative Mass n'étaient pas distribuées normalement. L'étape suivante est d'essayer de normaliser les données à l'aide de transformations mathématiques. Afin d'évaluer la normalité d'une variable, on trace un histogramme avec la fonction hist()
et on vérifie visuellement si la variable suit une distribution normale. Par exemple :
- | Vérifier la normalité des données avec la fonction hist()
# Grahpique Y ~ X avec une 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")
Une deuxième façon d'évaluer la normalité des données est d'utiliser le test de Shapiro-Wilk (fonction shapiro.test()
). Ce test compare la distribution des données observées à une distribution normale.
Les hypothèses nulle et contraire sont :
H0: les données observées sont distribuées normalement
H1: 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).
- Tester la normalité avec la fonction shapiro.test()
# 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
:
- Tester la normalité 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 indiquent tous que les variables ont besoin d'être transformées pour respecter la supposition de normalité (ex. une transformation logarithmique).
1.5 Transformation des données
Lorsque la supposition 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) | |
Asymétrie positive importante | log10(x) | |
Asymétrie positive importante | 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) où K est une constante soustraite de chaque valeur de x afin que la plus petite valeur soit 1 |
|
Asymétrie négative importante | log10(K - x) |
Dans notre cas, une transformation logarithmique (log10) 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é de nouveau.
- | Transformation de données
# 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) # Refaites 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)
1.6 Sortie du modèle
Lorsque les suppositions de base ont été vérifiées, les résultats du modèle peuvent être interprétés. On obtient ces résultats avec la fonction summary()
.
- | Sortie du modèle 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 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 associée apparraissent dans les deuxième et troisième colonnes respectivement. Donc,
β0 = 1.6724 ± 0.2472 est l'ordonnée à l'origine (± e.t.) du modèle de régression,
β1 = -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 valeurs de p associées sont les résultats d'un test statistique. Ce test vérifie si les coefficients de régression sont significativement différents de zéro. Dans notre cas, on peut voir que la variable logMass a une influence significative sur la variable logMaxAbund parce que la valeur de p associé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 la pente du modèle a une valeur négative.
Rappelez-vous qu'une corrélation entre deux variables n'implique pas nécessairement de relation de cause à effet. 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 R2 ajusté (ici 0.05484) et est calculé de la manière suivante :
où
p est le nombre total de paramètres de régression et n est la taille d'échantillon,
est la dispersion totale,
est la dispersion de la régression - aussi appelée la variance expliquée par le modèle.
Le R2 ajusté varie entre 0 et 1. Plus ce coefficient est élevé, meilleur est l'ajustment du modèle. 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 la valeur de p y étant 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.
1.7 Représentations graphiques
Les résultats d'une régression linéaire sont généralement représenté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 :
- | Tracer la régression Y ~ X avec une droite et des intervalles de confiance
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])
1.8 Sous-ensembles
Il est possible de réaliser une analyse sur seulement une partie des observations. Par exemple, on peut refaire l'analyse de régression en ne considérant que les oiseaux terrestres.
- | Régression sur un sous-ensemble d'observations
# 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 log10(MaxAbund) et log10(Mass) pour les passereaux (i.e. passerine birds).
Conseil :
La variable 'Passerine' est codée 0 et 1 comme la variable 'Aquatic'. Vous pouvez le vérifier avec la commande str(bird)
.
2. Test de t
Le test de t de Student, ou tout simplement test de t, permet de comparer les valeurs d'une variable réponse continue répartie entre deux groupes (ou traitements). Le test de t permet de déterminer si la moyenne d'un groupe est différente de celle de l'autre groupe. Le test de t est exprimé sous la forme suivante :
où
µ est la moyenne de la variable réponse,
correspond à l'effet du groupe i,
i prend la valeur de 1 ou 2,
correspond aux résidus du modèle (i.e. la variance inexpliquée).
Les hypothèses statistiques du test de t évaluent la différence entre les deux groupes :
H0: µ1 = µ2
H1: µ1 ≠ µ2
où
µ1 est la moyenne de la variable réponse pour le groupe 1,
µ2 est la moyenne de la variable réponse pour le groupe 2.
L'hypothèse contraire H1 affirme qu'il existe une différence entre les deux groupes au niveau de la variable réponse, ce qui constitue un test bilatéral.
Cependant, si le sens de la différence attendue est supportée par une hypothèse biologique, un test unilatéral peut être utilisé :
- si on s'attend à ce que la variable réponse soit supérieure pour le groupe 1, alors H1: µ1 > µ2,
- si on s'attend à ce que la variable réponse soit inférieure pour le groupe 1, alors H1: µ1 < µ2.
La statistique t du test de t qui est utilisée pour déterminer la valeur de p est calculée de la manière suivante :
où
1 et 2 sont les moyennes de la variable réponse y pour les groupes 1 et 2 respectivement,
s12 et s22 sont les variances de la variable réponse y pour les groupes 1 et 2 respectivement,
n1 et n2 sont les tailles d'échantillons des groupes 1 et 2 respectivement.
2.1 Suppositions de base
Si les suppositions de base du test de t ne sont pas respectées, les résultats du test peuvent être erronés. Ces suppositions concernent la forme de la distribution des données :
- 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 s12 et s22.
Pour ce faire, il faut effectuer un test de ratio des variances (i.e. un test de F).
Pour l'exemple ci-dessus, l'erreurs de type I est plus grande que la valeur α de l'échantillon #1. Il faut donc conclure qu'on ne peut pas rejeter l'hypothèse nulle, alors qu'en réalité, on aurait dû la rejeter !
Non-respect des suppositions
Si les variances entre les groupes ne sont pas égales, il est possible de corriger la situation avec la correction de Welch. Si les suppositions ne sont toujours pas respectées, il faut utiliser la version non paramétrique du test de t : 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 test de t apparié.
2.2 Effectuer un test de t
Dans R, les test de t sont exécutés avec la fonction t.test
. Par exemple, pour évaluer la différence de masse entre des oiseaux aquatiques et non aquatiques, vous devez utiliser le script suivant :
- | Test de t
# 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
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.
2.3 Effectuer un test de t avec la fonction lm()
Le test de t fait partie de la famille des modèles linéaires et est un cas spécifique de l'ANOVA (voir plus bas) à un critère et deux niveaux. Ceci signifie qu'on peut effectuer un test de t avec la fonction lm()
(qui signifie linear model) :
- | Le test de t comme modèle linéaire
ttest.lm1 <- lm(logMass ~ Aquatic, data=bird) anova(ttest.lm1) # La fonction anova permet de visualier les résultats du test
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
Lorsque les variances sont égales (i.e., test de t à deux échantillons), il est possible de démontrer que t2 = F:
- | Équivalence entre t et F
ttest1$statistic^2 anova(ttest.lm1)$F
2.4 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 doit écrire la commande de la façon suivante :
- | Test de t unilatéral
# Test de t en spécifiant l'argument "alternative" uni.ttest1 <- t.test(logMass~Aquatic, var.equal=TRUE, data=bird, alternative="less") 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.
3. ANOVA
L'analyse de la variance (ANOVA) est une généralisation du test de t de Student. L'ANOVA permet de voir si la moyenne d'une variable continue est différente entre trois ou plusieurs groupes ou traitements (Rappelez-vous que le nombre de groupes est limité à deux pour le test de t). L'ANOVA compare la variable réponse y entre les groupes (i.e. la variable explicative) de la manière suivante :
où
µ est la moyenne globale de la variable réponse,
Ai est l'effet du groupe i pour le facteur A,
i varie de 1 à n (n > 2),
εij sont les résidus du modèle (i.e. la variance inexpliquée).
L'ANOVA tente de détecter des différences au niveau de la variable réponse y entre les groupes en posant les hypothèses suivantes :
H0: µ1 = µ2 =… = µj =… = µn
H1: il y a au moins une moyenne µj différente des autres
L'ANOVA se base sur la partition de la somme des carrés des écarts à la moyenne pour déterminer si une hypothèse doit être acceptée ou rejetée. L'ANOVA 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.
La comparaison de la variance entre les traitements à celle intra-traitement permet de calculer la statistique F. Cette statistique correspond au ratio entre la moyenne des carrés des traitements (MSTrt) et la moyenne des carrés des erreurs (MSE). Ces deux termes sont obtenus en divisant leurs sommes des carrés respectives par leurs degrés de liberté (voir table ci-dessous). Une valeur de p peut ensuite être calculée à partir de la statistique de F qui suit une distribution de khi carré (χ2).
Source de variation | Degrés de liberté (df) | Somme des carrés des écarts à la moyenne | Moyenne des carrés | Statistique de F |
---|---|---|---|---|
Total | ||||
Facteur A | ||||
Résidus |
a: nombre de niveaux de la variable explicative A; r: nombre de répétitions par traitement; : moyenne globale de la variable réponse; i : moyenne de la variable réponse du traitement i.
3.1 Types d'ANOVA
- ANOVA à un critère de classification
Un facteur avec plus de deux niveaux - ANOVA à deux critères de classification (voir la section ci-dessous)
- Deux facteurs ou plus,
- Chaque facteur peut avoir de multiples niveaux,
- Les interactions entre chaque facteur 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 couvert dans cet atelier. Un modèle linéaire mixte peut également être utilisé pour ce type de données (voir l'atelier 6).
3.2 Suppositions de base
L'ANOVA doit respecter quelques suppositions statistiques pour que les résultats soient valides. Ces suppositions peuvent être vérifiées visuellement ou par des tests statistiques.
- Distribution normale
Les résidus d'un modèle d'ANOVA peuvent être visualisés à l'aide d'un diagramme quantile-quantile. 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 fonctionbartlett.test()
. Si la valeur de p de ce test est supérieure à 0.05, l'hypothèse nulle H0: s12 = s22 =… = sj2 =… = sn2 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.
Non-respect des suppositions
Si les suppositions ne sont pas respectées, vous pouvez essayer une transformation sur 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 transformer vos données, vous pouvez utiliser l'équivalent non paramétrique de l'ANOVA : le test de Kruskal-Wallis .
3.3 Contrastes
- 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. μ1 = μ2) ou des hypothèses plus complexes (e.g. (μ1 + μ2)/3 == μ3).
Le nombre de comparaisons doit être plus bas ou égal au nombre de degrés de liberté de l'ANOVA. Ces comparaisons doivent être indépendantes l'une de l'autre. Pour plus de détails, voyez la section avancée sur les contrastes plus bas.
3.4 Effectuer 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).
- | ANOVA
# 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)"))
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 :
3.5 Vérifications des suppositions
- Diagnostic du modèle
# 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.6 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.7 Tests complémentaires
Si l'hypothèse nulle est rejetée, il n'est pas possible de savoir quels niveaux de traitements sont différents des autres avec une simple ANOVA. Pour déterminer quels niveaux diffèrent des autres, il est nécessaire d'effectuer un test post hoc. Ce test compare les combinaisons de niveaux deux à deux et identifie où se trouvent les différences. 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 probablement celui qui est le plus utilisé. Dans R, on utilise la fonction TukeyHSD()
pour effectuer ce test :
- Test de Tukey
# À quel niveau se situe la différence de diète ? TukeyHSD(aov(anov1),ordered=T) # Cette comande est équivalente à la précédente : TukeyHSD(aov1,ordered=T)
La sortie de R retourne un tableau qui fait la liste de toutes les combinaisons deux à deux 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.8 Représentations graphiques
Après avoir vérifié les suppositions 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()
. Avec cette fonction, R produit un graphique de la variable réponse en fonction des niveaux de traitement. Les erreurs types ainsi que des lettres (représentant le résultat d'un test post hoc) peuvent y être ajoutées.
- Fonction barplot
# Représentation graphique 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)
3.9 Contrastes (section avancée et optionnelle)
4. ANOVA à deux critères de classification
Il est possible d'utiliser deux variables explicatives dans une ANOVA afin de mieux expliquer la variabilité d'une variable réponse (cf. section 3). Afin d'inclure une deuxième variable explicative dans une ANOVA, le modèle mathématique doit être réécrit de manière à inclure l'interaction entre ces deux variables :
où
µ est la moyenne globale de la variable réponse,
Ai est l'effet du niveau i du facteur A,
Bj est l'effet du niveau j du facteur B,
AiBj est l'interaction entre les deux facteurs,
i et j varient de 1 à n (n ≥ 2),
εijk sont les résidus du modèle.
Les hypothèses nulles d'une ANOVA à deux critères de classification sont légèrement différentes de celles d'une ANOVA à un critère de classification :
H01: Il n'y a pas de différence de moyenne parmi les niveaux du facteur A; µa1 = µa2 = … = µai =… = µan
H02: Il n'y a pas de différence de moyenne parmi les niveaux du facteur B; µb1 = µb2 = … = µbi =… = µbm
H03: Il n'y a pas d'interaction entre les facteurs A et B.
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 | ||||
Intra- cases (erreur) | ||||
Cases | ||||
Facteur A | ||||
Facteur B | ||||
Interaction entre A et B |
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
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.
4.2 Diagramme d'interaction
Les interactions peuvent être visualisées à l'aide de la fonction interaction.plot()
:
- Diagramme d'interaction
interaction.plot(bird$Diet, bird$Aquatic, bird$logMaxAbund, col="black", ylab = expression("log"[10]*"(Maximum Abundance)"), xlab="Diet")
Qu'est-ce que le trou sur la ligne des oiseaux aquatiques signifie ?
- Plan non équilibré
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 obtenir 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).
5. ANOVA non équilibrée (section avancée et optionnelle)
6. ANCOVA
L'analyse de covariance (ANCOVA) est un mélange de régression linéaire et d'ANOVA : on teste l'effet d'une (ou plusieurs) variable catégorique et d'une (ou plusieurs) variable continue sur une variable réponse continue. Le modèle mathématique sous-jacent de l'ANCOVA est définit de la manière suivante :
où
µ est la moyenne globale de la variable réponse,
Ai est l'effet du facteur A,
Bi est l'effet de la variable continue,
xij est la covariable mesurée sur l'observation yij,
i est la valeur moyenne de la covariable pour le groupe i,
i varie de 1 à n (n > 2) traitements,
εij sont les résidus du modèle.
Prenez note que le modèle est composé du terme Ai pour déterminer l'effet d'un traitement ou d'un facteur (comme dans une ANOVA) ainsi qu'un terme Βi pour tenir compte de l'effet d'une covariable (i.e. la pente d'une variable comme dans une régression linéaire). Donc, chaque traitement peut être décrit par une pente et une ordonnée à l'origine. En plus de tester si la variable réponse est influencée par au moins un niveau de la variable catégorique, l'ANCOVA teste également si la variable réponse est influencée par la variable continue (i.e. appelée covariable dans le cadre d'une ANCOVA). L'ANCOVA teste également si les niveaux de la variable catégorique influencent la variable réponse différemment en fonction de la valeur de la variable continue (i.e. l'interaction entre ces deux variables explicatives). Les hypothèses nulles de l'ANCOVA sont définies ainsi :
H01: Il n'y a pas d'effet de la variable catégorique (i.e. µ1 = µ2 =… = µi =… = µn)
H02: Il n'y a pas d'effet de la variable continue (i.e. β = 0)
H03: Il n'y a pas d'interaction entre la variable catégorique et la variable continue
6.1 Suppositions de base
Tout comme le test de t et l'ANOVA, l'ANCOVA doit respecter certaines suppositions statistiques qu'on peut vérifier à l'aide de diagrammes de diagnostic :
- Les résidus du modèle sont distribués normalement
- Homoscédasticité de la variance résiduelle
- Les résidus et les valeurs prédites sont indépendants,
- La variance résiduelle et les valeurs prédites sont indépendantes
- La variance égale entre les différents niveaux d'un facteur donné
- Les covariables ont toutes la même étendue de valeurs
- Les variables sont fixes
- Les facteurs et les covariables sont indépendants
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 modèles linéaires mixtes !
6.2 Types d'ANCOVA
Il est possible d'avoir plusieurs facteurs et variables au sein d'une même ANCOVA. Sachez cependant que 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 fréquentes 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-bas dans les deux derniers graphiques:
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 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 :
- Exemple d'ANCOVA
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 :
- Moyennes ajustées
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 afin de tester l'effet du facteur Diet, de la covariable Mass et de leur interaction sur la variable réponse MaxAbund.
7. Régression multiple
Une régression multiple teste les effets de plusieurs variables explicatives continues sur une variable réponse continue. La régression multiple se base sur le modèle mathématique suivant :
où
β0 est l'ordonnée à l'origine de la droite,
β1 est l'effet de la variable x1 (i.e. la pente de la droite de régression de la variable x1),
β2 est l'effet de la variable x2 (i.e. la pente de la droite de régression de la variable x2),
εi sont les résidus du modèle (i.e. la variance inexpliquée).
7.1 Suppositions de base
La régression multiple doit respecter certaines suppositions afin d'être valide :
- Distribution normale de la variable réponse
Ceci peut être vérifié par un test de Shapiro-Wilk (fonctionshapiro.test()
). Une transformation peut être utilisée si les données ne sont pas distribuées normalement.
- Orthogonalité
Les variables explicatives ne doivent pas être colinéaires (i.e. ne doivent pas être fortement corrélées). Si une variable explicative est corrélée avec une autre, ces deux variables vont probablement expliquer la même variabilité de la variable réponse : l'effet d'une variable va cacher l'effet de l'autre variable sur la variable réponse.
- Linéarité
Les relations entre les variables doivent être linéaires.
Les résidus du modèle doivent respecter les mêmes suppositions que celles de la régression linéaire simple, soient :
- Distribution normale des résidus
- Indépendance des résidus en fonction des variables explicatives xi (ou des valeurs prédites)
- Indépendance de la variance des résidus en fonction des variables explicatives xi (ou des valeurs prédites)
- Pas d'observations aberrantes (i.e. outliers)
Non-respect des suppositions
S'il existe une relation entre deux variables explicatives, elles sont 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 (voir l'atelier 9),
- Essayez une analyse pseudo-orthogonale.
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 en utilisant un modèle de régression multiple.
Défi 5
Est-il nécessaire d'appliquer une transformation à la variable réponse abund ?
Vous avez remarqué au défi 5 que la variable réponse abund n'a pas pu être normalisée. Ceci suggère qu'il faudrait peut-être laisser tomber la supposition de normalité et passer à un 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.
- Régression multiple
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")
7.3 Régression polynomiale (section avancée et optionnelle)
7.4 Régression pas à pas
Afin d'obtenir un modèle de régression multiple “optimal”, on peut débuter avec un modèle qui inclut toutes les variables explicatives et retirer les variables non significatives en procédant à 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 (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()
:
- Régression pas à pas
lm.full <- lm(abund ~ . - Present, data=Dickcissel) lm.step <- step(lm.full) summary(lm.step)
Les résultats indiquent que seulement six variables sont significatives parmi les treize 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.
7.5 Inflation de la variance
La colinéarité entre des variables explicatives peut être évaluée en calculant un facteur d'inflation de la variance à l'aide de la fonction vif()
de la librairie HH
:
- Facteur d'inflation de la variance
vif(clDD ~ clFD + clTmi + clTma + clP + grass, data=Dickcissel)
ce qui donne la sortie suivante :
clFD clTmi clTma clP grass 13.605855 9.566169 4.811837 3.196599 1.165775
On considère qu'un facteur d'inflation de la variance supérieur à cinq indique que les variables sont colinéaires. La sortie de R indique que les variables clDD, clFD et clTmi sont fortement colinéaires. Seulement une des trois variables doit être retenue dans le modèle de régression final.
8. Partition de la variation (section optionnelle et avancée)
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 suppositions 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