Loading presentation...

Present Remotely

Send the link below via email or IM

Copy

Present to your audience

Start remote presentation

  • Invited audience members will follow you as you navigate and present
  • People invited to a presentation do not need a Prezi account
  • This link expires 10 minutes after you close the presentation
  • A maximum of 30 users can follow your presentation
  • Learn more about this feature in our knowledge base article

Do you really want to delete this prezi?

Neither you, nor the coeditors you shared it with will be able to recover it again.

DeleteCancel

CSBQ Atelier R 4

No description
by

CSBQ QCBS

on 19 October 2016

Comments (0)

Please log in to add your comment.

Report abuse

Transcript of CSBQ Atelier R 4

Centre de la science de la biodiversité du Québec

Série d'ateliers R

Atelier 4:
Modèles linéaires

Site web
:
http://qcbs.ca/wiki/r/atelier4

Objectif d'enseignement
Télécharger les données
birdsdiet
:
setwd()
bird<-read.csv("birdsdiet.csv")
Type de variable explicative
1 variable
1. Régression
linéaire simple

2 ou plus
5. Régression
linéaire multiple

1 avec 2 niveaux
2. test t
3. ANOVA
Quantitative
Qualitative
Quantitative ET qualitative
4. ANCOVA
Nous voulons d'abord vérifier si l'abondance maximale des oiseaux (MaxAbund) est une fonction de la masse des oiseaux (Mass)
Variable réponse:
abondance d'oiseaux /
num: quantitative
Variable explicative :
masse /
num : quantitative
Dans R, la fonction
lm()
est utilisée pour ajuster un modèle linéaire
lm1 <- lm(Y ~ X)
Variable réponse
Nom de l'objet contenant le modèle linéaire
Exécution du modèle linéaire dans R
Suppositions
Régression linéaire simple

Yi = α + βXi + εi

Relation linéaire entre la variable réponse (Y) et explicative (X)
Variable explicative continue
Distribution normale et homogène des résidus
Indépendance des résidus
Aucune observation aberrante
Étape 1:
Créer votre modèle linéaire
Étape 2:
Vérifier les suppositions
Suppositionsrespectées
Étape 3:
Estimer les paramètres, le seuil de signification des paramètres, tracer le modèle
Suppositions non-respectées
Pouvez-vous transformer vos variables (est-ce justifié)?
Retourner à l'
étape 1
avec des variables transformées
Oui
Non
Essayer GLM qui pourraient mieux convenir aux données
Étape 1:
Créer votre modèle linéaire
par()
définit les paramètres du graphique, par exemple, l'argument
mfrow
spécifie le nombre de rangées et colonnes
plot()
est la fonction pour faire le graphique
par(mfrow=c(2,2))
plot(lm1)
La sortie comprend les quatre graphiques diagnostics de la fonction lm()
Devrait montrer une dispersion de points sans patron
Exemple d'indépendance
(ce que nous recherchons!)
Exemples de non-indépendance

(ce que nous ne voulons pas!)
Solution :
Transformer vos données ou essayer une distribution autre que linéaire (gaussienne); modèle linéaire généralisé (GLM) (ex: Poisson, binomial, binomial négatif, etc.)
Graphique # 1 - Résidus vs valeur ajustée
Si les points se situent de façon linéaire sur la ligne 1: 1, les résidus suivent une distribution normale
Graphique # 3 - Normal QQ
Graphique # 2 - échelle localité
Devrait montrer une dispersion de points sans patron
Aucun patron dans les résidus
Forte tendance dans les résidus
Compare la distribution (quantiles) des résidus aux quantiles d'une distribution normale
Recherche les valeurs influentes
Points de levier :
observations extrêmes ou périphériques de la variable explicative. Parce qu'ils n'ont pas d'observations voisines, la ligne de régression passe près de ces points.
Ils peuvent (ou pas) avoir une grande influence sur la régression
Graphique # 4 - Résidus vs effet de levier
Effet levier vs influence
Les points de levier avec une forte influence peuvent être identifiés avec une
distance de Cook supérieure à 0,5

Pas d'effet de levier
Faible influence
Effet de levier
Pas d'influence
Effet de levier
Influence élevée
Aucune valeur influente
Effet de levier élevé et influence raisonnable
Ici, le point 32 a un effet de levier

élevé, mais son influence est acceptable (à l'intérieur des limites de la distance Cook de 0,5)
Effet de levier et influence élevée
Points en dehors de la limite de 0,5 de la distance Cook
Ces points ont trop d'influence sur la régression
La fonction
summary()
est utilisée pour obtenir les paramètres, leur importance, etc :
Étape 3 :
Estimer les paramètres et leur seuil de signification
** Vous ne devriez jamais supprimer les valeurs aberrantes si vous n'avez pas de bonnes raisons de le faire (ex: erreur de mesure)
GLM
GLM
Discussion de groupe
Défi 1
Examiner la relation entre logMaxAbund et logMass chez les passereaux ("passerine birds")

Sauvegarder l'objet du modèle sous lm4
INDICE: comme les espèces aquatiques, les passereaux sont codées 0/1, ce qui peut être vérifié à partir de la structure de la base de données
Rappel : test de t
Variable réponse
: quantitative
Variable explicative
: qualitative avec 2 niveaux
Suppositions
Les données suivent une distribution normale
Les variances des groupes sont homogènes

Test de Mann-Whitney

-

L'équivalent
non paramétrique
du test de t lorsque les suppositions ne sont pas respectées
Exécuter un test de t dans R
t.test(Y~X2, data= data, alternative = "two.sided")
Vous pouvez utiliser la fonction
t.test()
variable
réponse
facteur
(2 niveaux)
nom du jeu de données
hypothèse alternative :
"two.sided" (par défaut), "less", ou "greater"
Le test de t est un modèle linéaire et un cas spécifique de l'ANOVA avec un facteur à 2 niveaux
lm.t <-lm(Y~X2, data = data)
anova(lm.t)
moyennes des 2 groupes
Exécuter un test t dans R
Les oiseaux aquatiques sont-ils plus lourds que les oiseaux terrestres?
Variable réponse :
Bird mass
/

num: quantitative
Variable explicative :
Aquatic
/

2 niveaux: 1 or 0 (oui or non)
Exécuter un test de t dans R
1 avec 3 niveaux ou plus
et/ou plusieurs facteurs

Analyse de Variance (ANOVA)
Décomposition de la variation observée de la variable réponse en effets additifs d'un ou de plusieurs facteurs et de leurs interactions
Y = μ +
+ Interactions between factors

+ Residuals

moyenne globale de la variable réponse
Rappel : ANOVA
Suppositions
Normalité des résidus
L'égalité de la variance inter-groupes
Test complémentaire
Lorsque l'ANOVA détecte une différence significative entre les groupes, l'analyse n'indique pas quel(s) groupe(s) diffère(nt) de(s) l'autre(s)
Un test couramment utilisé
a posteriori
pour répondre à cette question est le

Test de Tukey
Rappel : ANOVA
Comme le test de t, l'ANOVA est plus robuste lorsque la taille de l'échantillon est grande et les groupes ont des tailles égales
Pour chaque facteur et chaque terme d'interaction :
H0 : La moyenne de la variable réponse est égale entre les groupes
H1 : Au moins deux groupes ont une moyenne différente
Compare la variation intra-groupe et inter-groupe afin de déterminer si les moyennes des groupes diffèrent
Généralisation du test t à >2 groupes, et / ou ≥ 2 facteurs explicatifs
Exécuter une ANOVA dans R
Il est de nouveau possible d'utiliser la fonction

lm()
anov1 <- lm(logMaxAbund ~ Diet, data=bird)
anova(anov1)
Il y a aussi une fonction spécifique pour l'analyse de la variance dans R,

aov()
aov1 <- aov(logMaxAbund ~ Diet, data=bird)
summary(aov1)
Quand un terme d'interaction est significatif, cela signifie que l'effet d'un facteur varie en fonction des niveaux d'un autre facteur
Vérifier les suppositions
Test de Bartlett

-

égalité de la variance entre les groupes
bartlett.test(logMaxAbund ~ Diet, data=bird)
Test de Shapiro-Wilk

- Normalité des résidus
shapiro.test(resid(anov1))
Et si les suppositons ne sont pas respectées...
Transformer vos données
- pourrait égaliser les variances et normaliser les résidus, et peut convertir un effet multiplicatif en un effet additif
data$logY <- log10(data$Y)
* ré-exécuter votre modèle avec la variable transformée et vérifier à nouveau les hypothèses
Test de
Kruskal-Wallis

- équivalent non paramétrique de l'ANOVA si vous ne pouvez pas
(ou ne voulez pas)
transformer vos données
kruskal.test(Y~X, data)
* voir le wiki de l'atelier 1 pour les règles de transformation de données
$ Family : Factor w/ 53 levels "Anhingas","Auks& Puffins",..: 18 25 23 21 2 10 1 44 24 19 ...
$ MaxAbund : num 2.99 37.8 241.4 4.4 4.53 ...
$ AvgAbund : num 0.674 4.04 23.105 0.595 2.963 ...
$ Mass : num 716 5.3 35.8 119.4 315.5 ...
$ Diet : Factor w/ 5 levels "Insect","InsectVert",..: 5 1 4 5 2 4 5 1 1 5 ...
$ Passerine: int 0 1 1 0 0 0 0 0 0 0 ...
$ Aquatic : int 0 0 0 0 1 1 1 0 1 1 ...
Est-ce que l'abondance maximale dépend du régime alimentaire ?
Variable réponse :
MaxAbund
/

num : quantitative
Variable explicative :
Diet
/

facteur avec 5 niveaux
Exécuter une ANOVA dans R
Sorties de notre modèle ANOVA
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
Triage en ordre alphabétique des niveaux et comparaison au niveau de référence ("Insect")
Différence significative entre les groupes, mais nous ne savons pas lesquels!
Test a posteriori
TukeyHSD(aov(anov1),ordered=T)
$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

Seuls "Vertebrate" et "PlantInsect" diffèrent
Interprétation des contrastes
Défi 2
Testez si l'abondance maximale varie à la fois en fonction du régime alimentaire (Diet) et de l'habitat (Aquatic)
Test de t apparié
- Lorsque

les deux groupes ne sont
pas indépendants
(par exemple, des mesures sur la même personne récoltées lors de 2 années différentes)
Rappel : Test t
Correction de Welch

- Lorsque les écarts entre les groupes ne sont pas égaux

(option par défaut dans R!)
Types d'ANOVA
ANOVA à un critère de classification
Un facteur avec plus de 2 niveaux
ANOVA à deux critères de classification
2 facteurs ou plus
Chaque facteur peut avoir plusieurs niveaux
Vous devez tester les
interactions
Les modèles linéaires à effets mixtes peuvent également être utilisés pour ce type de données (atelier LMM)
Mesures répétées?
ANOVA à deux critères de classification
anova.avec.lm <- lm(Y~X*Z, data)
lorsque vous utilisez le symbole "*" avec lm(), le modèle inclut les effets de chaque facteur séparément,
ainsi que leur interaction
Analysis of Variance Table

Response: Y
Df Sum Sq Mean Sq F value Pr(>F)
X 4 5.1059 1.27647 3.0378 0.02669 *
Z 1 0.3183 0.31834 0.7576 0.38870
X:Z 3 2.8250 0.94167 2.2410 0.10689
Residuals 45 18.9087 0.42019
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Exemple d'interaction non-significative
anova.avec.lm <- lm(Y~X+Z, data)
lorsque vous utilisez le symbole "+" avec lm(), le modèle inclut les effets de
chaque facteur séparément
(pas d'interaction)
Selon le principe de
parcimonie
, vous voulez que votre modèle explique le plus possible de la variance observée dans les données, avec le moins de termes possible

Enlever le terme d'interaction s'il n'est pas significatif, et ré-exécuter le modèle

Analyse de covariance (ANCOVA)
Sous-ensembles de la variation observée :
Y = μ +
+ Interactions entre facteurs
Effets principaux
des facteurs
+ Résidus

You will see in the assumptions that there should be NO interaction between covariates and factors
ANCOVA avec une covariable et un facteur
1- Déterminer l'effet du facteur et de la covariable sur ​​la variable réponse
2- Déterminer l'effet du facteur sur la variable réponse après avoir enlevé l'effet de la covariable
3- Déterminer l'effet de la covariable sur la variable réponse en contrôlant l'effet du facteur
Objectifs de l'analyse
Suppositions
Combinaison de l'ANOVA et de la régression linéaire
ancova.exemple <- lm(Y~X*Z, data=data)
summary(ancova.exemple)
Variables sont
fixes*

Distribution normale des résidus
Indépendance des résidus
Variance égale entre les différents niveaux d'un facteur
Aucune valeur aberrante
Aucune interaction entre les facteurs et covariable(s)
Types d'ANCOVA
ANCOVA fréquemment utilisées
Vous pouvez avoir n'importe quel nombre de facteurs et / ou variables, mais lorsque leur nombre augmente, l'interprétation des résultats devient de plus en plus complexe
1- Une
covariable et un facteur
2- Une covariable et deux facteurs
3- Deux covariables et un facteur
Les variables explicatives sont un mélange de variables quantitatives (covariable) et qualitatives (facteurs)
+ Effets principaux
des covariables
+ Interactions entre covariables et facteurs
Si non-respectées, vous pouvez essayer de transformer vos données
2- Vérifiez si votre interaction est significative
ancova.exemple2 <- lm(Y~X+Z, data=data)
Si vous avez une interaction significative entre votre facteur et votre covariable, vous ne pouvez pas atteindre ces objectifs!
Pas d' interaction
Interaction :
Un niveau du facteur a une pente différente
Interaction :
De nombreux niveaux ont des pentes différentes
1- Exécutez un modèle pour tester les effets du régime alimentaire (Diet), de la masse (logMass) ainsi que leur interaction sur l'abondance maximale des oiseaux (logMaxAbund)
Si votre covariable et votre facteur sont significatifs, vous avez un cas comme celui-ci:
summary(ancova.exemple2)
Si vous voulez comparer les moyennes des différents facteurs, vous pouvez utiliser les
moyennes ajustées

La fonction

effect()
utilise les équations données par l'ANCOVA pour estimer les moyennes de chaque niveau, corrigées pour l'effet de la covariable
ancova.exemple <- lm(Y ~ X*Z, data=data)
# X = variable quantitative; Z = variable qualitative

library(effects)
adj.means.ex <- effect('Z', ancova.exemple)
plot(adj.means.ex)
* Une variable
fixe
est celle à laquelle vous êtes particulièrement intéressés (c.-à-d. masse de l'oiseau). En revanche, une variable
aléatoire
représente du bruit que nous souhaitons contrôler (c.-à-d. le site où un oiseau a été échantillonné)
> Voir l'atelier sur les LMMs
Si seulement votre facteur est significatif, éliminer la covariable -> vous avez une
ANOVA
Si votre interaction covariable * facteur est significative, vous voudrez peut-être tester quel(s) niveau(x) du facteur a(ont) des pentes différentes.
Si seulement votre covariable est significative, éliminer le facteur -> vous avez une
régression linéaire simple
Vérifier vos suppositions !
Très similaire à ce que vous avez fait précédemment
Nous ne considérerons que le premier cas aujourd'hui, mais les deux autres sont similaires
$ Family : Factor w/ 53 levels "Anhingas","Auks& Puffins",..: 18 25 23 21 2 10 1 44 24 19 ...
$ MaxAbund : num 2.99 37.8 241.4 4.4 4.53 ...
$ AvgAbund : num 0.674 4.04 23.105 0.595 2.963 ...
$ Mass : num 716 5.3 35.8 119.4 315.5 ...
$ Diet : Factor w/ 5 levels "Insect","InsectVert",..: 5 1 4 5 2 4 5 1 1 5 ...
$ Passerine: int 0 1 1 0 0 0 0 0 0 0 ...
$ Aquatic : int 0 0 0 0 1 1 1 0 1 1 ...
L'abondance maximale varie-t-elle en fonction du régime alimentaire et la masse des oiseaux?
Variable réponse:
MaxAbund
/

num : quantitative continue
Variables explicatives:
Diet
/

facteur à 5 niveaux
Mass
/ numérique continue
Rappel: régression linéaire multiple
Variables explicatives
: 2 ou plus variables continues
Il s'agit d'une généralisation de la régression linéaire simple
Suppositions
Les variables explicatives suivent une distribution normale
La relation en chaque variable explicative et la variable réponse est de type linéaire
Les variables explicatives sont orthogonales (pas de colinéarité)
Normalité des résidus
Indépendance des résidus
Pas de valeur aberrante

Yi = α + β1X1i + β2X2i + β3X3i + ...+ βkXki+ ε
*Comme la régression linéaire simple!
La colinéarité

Vérifiez les suppositions
lm.mult <- lm(abund ~ clTma + NDVI + grass, data=Dickcissel)
summary(lm.mult)
plot(data=Dickcissel)
Si vous observez un patron entre vos deux variables explicatives, elles peuvent être colinéaires!
Vous devez éviter ceci, sinon leurs effets sur la variable réponse seront confondus



Gardez une seule des variables colinéaires
Essayez l'analyse multivariée (
ateliers 9 et 10
)

Vérifiez les autres suppositions, comme pour la régression linéaire simple

par(mfrow=c(2,2))
plot(lm.mult)
Souvenez-vous du principe de parcimonie: expliquer le plus de variation avec le plus petit nombre de termes dans votre modèle

Quel est le meilleur modèle?
En utilisant le jeu de données
Dickcissel
,
comparez l'importance relative du climat (clTma), de la productivité (NDVI) et de la couverture du sol (grass) comme prédicteurs de l'abondance de dickcissels (abund)
Étape 2:
Vérifier les suppositions avec les graphiques diagnostics
lm1 <- lm(bird$MaxAbund ~ bird$Mass)
bird$logMaxAbund <- log10(bird$MaxAbund)
bird$logMass <- log10(bird$Mass)
Comparer la variance expliquée par lm2, lm3 et lm4
Pouvons-nous améliorer le modèle si nous analysons que les oiseaux terrestres?
lm3 <- lm(logMaxAbund~logMass, data=bird, subset=!bird$Aquatic)
# élimine les oiseaux aquatiques (= TRUE)

# ou également
lm3 <- lm(logMaxAbund~logMass, data=bird, subset=bird$Aquatic == 0)

# Examiner les graphiques diagnostic
par(mfrow=c(2,2))
plot(lm3)
summary(lm3)

# Comparer les deux modèles de régression
par(mfrow=c(1,2))
plot(logMaxAbund~logMass, data=bird)
plot(logMaxAbund~logMass, data=bird, subset=!bird$Aquatic)
ttest1 <- t.test(logMass~Aquatic, var.equal=TRUE, data=bird)

# Ou utilisez lm()
ttest.lm1 <- lm(logMass ~ Aquatic, data=bird)
boxplot(logMaxAbund ~ Diet, data=bird)
INDICE: Examinez les facteurs Diet, Aquatic et leur interaction avec une ANOVA à deux critères de classification

lm(Y ~ A*B)
où A est le premier facteur, B le deuxième et "*" décrit l'interaction
opar <- par(mfrow=c(1,3))
plot(abund ~ clTma, data=Dickcissel)
plot(abund ~ NDVI, data=Dickcissel)
plot(abund ~ grass, data=Dickcissel)
par(opar)
Comme nous l'avons remarqué dans la section sur la
régression linéaire multiple
, certaines variables semblent avoir des relations non-linéaires avec la variable MaxAbund.
Pour tester des relations non-linéaires, des régressions polynomiales de différents degrés sont comparées.
Régression polynomiale
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)
Régression polynomiale
Comparez les modèles polynômiaux et déterminez quel modèle niché nous devrions sélectionner.
Exécutez un résumé de ce modèle, reportez l'équation de la régression, les valeurs de p, et le R carré ajusté

Régression pas à pas
Pour choisir le meilleur modèle, nous allons soustraire successivement un terme différent du modèle complet et comparer tous les modèles entre eux
La fonction
step()
soustrait un terme au modèle de façon itérative et sélectionne le meilleur modèle
c-à.d. le modèle avec le Critère d'Information Akaike (AIC) le plus bas
# Régression multiple contenant toutes les variables sauf la variable "Present"
lm.full <- lm(abund ~ . - Present, data=Dickcissel)
summary(lm.full)

lm.step <- step(lm.full)
summary(lm.step)
?step
Partitionnement de la variation
Par défaut, "contr.treatment" compare chaque niveau du facteur à un niveau de référence.
L'estimation de l'ordonnée à l'origine (ici = 1.1539) est le niveau de référence et correspond à la moyenne du premier niveau (en ordre alphabétique) du facteur "Diet".
# Calculez l'ordonnée à l'origine de référence + l'ordonnée à l'origine de chaque niveau de Diet


tapply(bird$logMaxAbund, bird$Diet, mean)
coef(anov1)
coef(anov1)[1] + coef(anov1)[2] # InsectVert
coef(anov1)[1] + coef(anov1)[3] # Plant
Que remarquez-vous?
Interprétation des contrastes
# Comparez le niveau "Plant" à tous les autres niveaux du facteur Diet
bird$Diet2 <- relevel(bird$Diet, ref="Plant")
anov2 <- lm(logMaxAbund ~ Diet2, data=bird)
summary(anov2)
anova(anov2)


# Ordonner les niveaux selon leur médiane
bird$Diet2 <- factor(bird$Diet, levels=names(med))
anov2 <- lm(logMaxAbund ~ Diet2, data=bird)
summary(anov2)
anova(anov2)
Il se peut que vous vouliez définir un niveau de référence différent
Observez-vous un changement quant aux niveaux du facteur "Diet" qui sont significatifs?
Interprétation des contrastes
Un point important à remarquer à propos du contraste par défaut dans R ("contr.treatment") est qu'il n'est PAS orthogonal.
# Pour être orthogonal, les propriétés suivantes doivent être respectées:
# 1. La somme des coefficients égale 0
sum(contrasts(bird$Diet)[,1])
# 2. La somme du produit de deux colonnes égale 0
sum(contrasts(bird$Diet)[,1]*contrasts(bird$Diet)[,2])
# La propriété 1 n'est pas respectée

# Changez les contrastes pour mettre les niveaux orthogonaux
options(contrasts=c("contr.helmert", "contr.poly"))
sum(contrasts(bird$Diet)[,1])
# Propriété 1
sum(contrasts(bird$Diet)[,1]*contrasts(bird$Diet)[,2])
# Propiété 2
anov3 <- lm(logMaxAbund ~ Diet, data=bird)
summary(anov3)
# Les contrastes Helmert vont contraster le deuxième niveau avec le premier, le troisième avec la moyenne des deux premiers niveaux, etc.

Certaines variables explicatives de la
régression linéaire multiple
étaient fortement corrélées(c.-à-d.
multicolinéarité
)
Ces variables devraient être enlevées du modèle.
La colinéarité entre variables explicatives peut être détectée à l'aide de critères d'inflation de la variance (fonction vif du progiciel "HH").
vif(clDD ~ clFD + clTmi + clTma + clP + grass, data=Dickcissel)
# les valeurs supérieures à 5 sont considérées colinéaires
.

Partitionnement de la variation
Utilisez varpart() afin de partitionner la variation de la variable
abund
avec toutes les variables de la couverture du paysage groupées ensemble et toutes les variables du climat groupées ensemble (laissez NDVI à part).
Note
: les variables colinéaires n'ont pas besoin d'être enlevées avant l'analyse.
library(vegan)
part.lm = varpart(Dickcissel$abund,
Dickcissel[ ,c("clDD","clFD","clTmi","clTma","clP")],
Dickcissel[ ,c("broadleaf","conif","grass","crop", "urban","wetland")])
part.lm

showvarparts(2)
plot(part.lm,digits=2,bg=rgb(48,225,210,80,maxColorValue=225),col="turquoise4")
Partitionnement de la variation
la proportion de la variation de la variable abund expliquée par le climat seulement est 28.5% (obtenu par X1|X2), par la couverture du paysage seulement est ~0% (X2|X1), et par les deux combinés est 2.4%.
# Tester si chaque fraction est significative

# Climat seul
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)

# Couverture du paysage seul
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)

# Conclusion: la fraction expliquée par la couverture du paysage n'est pas significative une fois que nous avons pris en compte l'effet du climat.

Partitionnement de la variation
R2-adj a changé de 0,05 à 0,25 quand nous avons exclu les oiseaux aquatiques :
Exécuter un test de t dans R
boxplot(logMass ~ Aquatic, data=bird,names=c("Non-Aquatic","Aquatic"))
Premièrement, visualisation les données à l'aide de la fonction
boxplot()
Représentation graphique de l'ANOVA à l'aide de la fonction
barplot()
Représentation graphique
Il est important de noter que la variable réponse ne varie pas de façon linéaire avec les variables explicatives
Section avancée
1. Interprétation des contrastes
2. ANOVA non équilibrée
table(bird$Aquatic)

unb.anov1 <- lm(logMaxAbund ~ Aquatic + Diet, data=bird)
unb.anov2 <- lm(logMaxAbund ~ Diet + Aquatic, data=bird)
anova(unb.anov1)
anova(unb.anov2)
# L'ordre des covariables a changé la valeur de la somme des carrés.


# Maintenant essayez une ANOVA de type III
library(car)
Anova(unb.anov1,type="III")
Anova(unb.anov2,type="III")

ANOVA non-équilibrée
Le jeu de données Birdsdiet est en réalité non équilibré (le nombre d'espèces aquatiques n'égale pas le nombre d'espèces non-aquatiques).
Que remarquez-vous en utilisant Anova()?
4. Partitionnement de la variation
Hypothèses
H0
: Moyenne du groupe 1 est égale à moyenne du groupe 2
H1
: Moyennes des deux groupes diffèrent
Si p<0,01 (ou 0,05), l'hypothèse de l'absence de différence entre les moyenne des 2 groupe (H0) peut être rejetée, avec un risque de 0,01 (ou 0,05) de se tromper.
Exécuter un test de t dans R
Testons l'homogénéité des variances avec la fonction
var.test()
var.test(logMass~Aquatic,data=bird)
Le rapport des variances n'est pas statistiquement différent de 1, celles-ci peuvent donc être considérées comme égales
Discussion de groupe
Les oiseaux aquatiques sont-ils plus lourds que les oiseaux terrestres ?
# Test t unilatéral
uni.ttest1 <- t.test(logMass~Aquatic, var.equal=TRUE, data=bird, alternative="less")
uni.ttest1
Exécuter un test de t dans R

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 in group 0 mean in group 1
1.583437 2.908289
Il existe une différence entre la masse des oiseaux aquatiques et terrestres
Qu'avez-vous conclu ?
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
Oui, les oiseaux aquatiques sont plus lourds que les oiseaux terrestres
Discussion de groupe
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)
segments(bp, means - se, bp, means + se, lwd=2)
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)
Les deux tests sont non-significatifs ; les résidus du modèle peuvent être considérés normaux et les variances homogènes
* Voir la
section avancée
sur la
régression polynomiale
pour la solution!
3. Régression polynomiale
Optionnel
En utilisant le jeu de données Dickcissel, testez la relation non-linéaire entre l'abondance et la température en comparant trois modèles polynômiaux groupés (de degrés 0, 1, and 3):
Un modèle polynômial ressemble à ceci:
termes
Ce polynôme a trois termes
Régression polynomiale
Pour un polynôme avec une variable (comme x), le degré est l'exposant le plus élevé de cette variable.
Nous avons ici un polynôme de degré 4
Régression polynomiale
Lorsque vous connaissez le degré, vous pouvez lui donner un nom:
Régression polynomiale
*Voir le script R avancé (optionnel)pour la solution
Traçons le graphique Y ~ X avec la droite de régression et les histogrammes de Y et X pour explorer leurs distributions
Normalisons les données en appliquant une transformation log10()
Ajoutons ces variables transformées à notre base de données
Vérifions la normalité des données à l'aide d'un test de
Shapiro-Wilk
et d'un test d'asymétrie (
skewness
)
shapiro.test(bird$MaxAbund)
shapiro.test(bird$Mass)
skewness(bird$MaxAbund)
skewness(bird$Mass)
Shapiro-Wilk normality test
data: bird$MaxAbund
W = 0.5831, p-value = 3.872e-11

Shapiro-Wilk normality test
data: bird$Mass
W = 0.5467, p-value = 1.155e-11
[1] 3.167159
# mesure de symétrie de MaxAbund
[1] 3.146062
# mesure de symétrie de Mass
Dans les deux cas, les distributions ne sont pas normales
La valeur positive indique que la distribution des données est décalée vers la gauche
Beaucoup mieux!
Pouvez-vous écrire l'équation de la droite de régression pour votre modèle lm2?

Les paramètres sont-ils importants?

Quelle proportion de la variance est expliquée par le modèle lm2?
plot(lm2)
Défi 1 - Solution
# Exécutons la régression linéaire
lm4 <- lm(logMaxAbund ~ logMass, data=bird, subset=bird$Passerine == 1)

# Examiner les graphiques diagnostics
par(mfrow=c(2,2))
plot(lm4)
summary(lm4)

# Comparer la variance expliquée par lm2, lm3 et lm4
str(summary(lm4))
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

Le meilleur modèle parmi les trois est lm3
(seulement les oiseaux terrestres)
Call:
lm(formula = abund ~ clDD + clFD + clTmi + clTma + clP + grass,
data = Dickcissel)

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
Variables sélectionnées par step()
(... si le temps le permet)
Vérifiez que
t.test()

et

lm()
donnent le même modèle
ttest1$statistic^2
anova(ttest.lm1)$F
# réponse : F=60.3845 dans les deux cas
Les 3 variables sont importantes. On garde tout!
Call:
lm(formula = abund ~ clTma + NDVI + grass, data = Dickcissel)

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
Nous pouvons maintenant procéder au test t!
Lorsque la supposition d'égalité de variance est confirmée, t^2 = F
Effets principaux
des facteurs
Y = u +
+ Résidus
+ Interactions
entre facteurs
Visualisons tout d'abord les données avec la fonction
boxplot()

(par défaut en ordre alphabétique)
Nous pouvons changer l'ordre des niveaux afin qu'il suivent l'ordre croissant de leurs médianes respectives en utilisant les fonctions
tapply()
et
sort()
med <- sort(tapply(bird$logMaxAbund, bird$Diet, median))
boxplot(logMaxAbund ~ factor(Diet, levels=names(med)))
Enlevez la variable qui est la moins significative

Exécution du modèle linéaire dans R
$ Family : Factor w/ 53 levels "Anhingas","Auks& Puffins",..: 18 25 23 21 2 10 1 44 24 19 ...
$ MaxAbund : num 2.99 37.8 241.4 4.4 4.53 ...
$ AvgAbund : num 0.674 4.04 23.105 0.595 2.963 ...
$ Mass : num 716 5.3 35.8 119.4 315.5 ...
$ Diet : Factor w/ 5 levels "Insect","InsectVert",..: 5 1 4 5 2 4 5 1 1 5 ...
$ Passerine: int 0 1 1 0 0 0 0 0 0 0 ...
$ Aquatic : int 0 0 0 0 1 1 1 0 1 1 ...
plot(lm1)
summary(lm2)
Coefficients:
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
lm2$coef
str(summary(lm2))
summary(lm2)$coefficients
summary(lm2)$r.squared
Nous pouvons aussi extraire les paramètres du modèle, par exemple :
Graphique # 1 - Résidus vs valeur ajustée
Régression linéaire multiple dans R
Dickcissel <- read.csv("Dickcissel.csv")
str(Dickcissel)
Visualisez le tableau de la structure des données en utilisant la fonction
str()
:
str(bird)
Exécution du modèle linéaire dans R
$ Family : Factor w/ 53 levels "Anhingas","Auks& Puffins",..: 18 25 23 21 2 10 1 44 24 19 ...
$ MaxAbund : num 2.99 37.8 241.4 4.4 4.53 ...
$ AvgAbund : num 0.674 4.04 23.105 0.595 2.963 ...
$ Mass : num 716 5.3 35.8 119.4 315.5 ...
$ Diet : Factor w/ 5 levels "Insect","InsectVert",..: 5 1 4 5 2 4 5 1 1 5 ...
$ Passerine: int 0 1 1 0 0 0 0 0 0 0 ...
$ Aquatic : int 0 0 0 0 1 1 1 0 1 1 ...
Étape 2:
Vérifier les suppositions de
lm1
Transformer les données
plot(logMaxAbund ~ logMass, data=bird)
abline(lm2)
hist(log10(bird$MaxAbund))
hist(log10(bird$Mass))
Étape 1:
Exécuter une régression linéaire sur les données transformées
lm2 <- lm(bird$logMaxAbund ~ bird$logMass)
Étape 2:
Vérifier les suppositions de notre modèle
lm2
Étape 2:
Vérifier les suppositions du modèle
lm2
plot(bird$MaxAbund ~ bird$Mass)
abline(lm1)
hist(bird$MaxAbund)
hist(bird$Mass)
* Le test est plus robuste lorsque la taille de l'échantillon est plus élevée et lorsque les groupes ont des tailles égales
Vous pouvez donc aussi utiliser la fonction
lm()
Exécuter une ANOVA dans R
summary.lm(aov1)
summary(anov1)
Lorque l'ANOVA détecte un effet significatif de la variable explicative, un test post-hoc, avec la fonction
TukeyHSD()
,
doit être effectué pour déterminer quel(s) traitement(s) diffère(nt)

ANOVA à deux critères de classification
Défi 2 - Solution
anov2 <- lm(logMaxAbund ~ Diet*Aquatic, data=bird)
summary(anov2)
anova(anov2)
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
ANCOVA avec une covariable et un facteur
ANCOVA avec une covariable et un facteur
Défi 3
ANCOVA avec une covariable et un facteur
ANCOVA avec une covariable et un facteur
Variable indépendante
Spécifier le répertoire (fichier)
graphique 2-par-2
Suppositions non respectées - quelle est la cause?
Discussion de groupe
**Vous pouvez exclure des objets en utilisant "=!"
F test to compare two variances
data: logMass by Aquatic
F = 1.0725, num df = 38, denom df = 14, p-value = 0.9305
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.3996428 2.3941032
sample estimates:
ratio of variances
1.072452
Spécifie que l'homogénéité des variances est respectée
uni.ttest1
Défi 3: Solution
ancov1 <- lm(logMaxAbund ~ logMass*Diet, data=bird)
anova(ancov1)
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
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Interaction entre logMass et Diet est non-significative
ancov2 <- lm(logMaxAbund ~ logMass + Diet, data=bird)
anova(ancov2)
Éliminer le terme d'interaction, puis ré-évaluer le modèle contenant les effets simples de logMass et Diet
Défi 3: Solution
Analysis of Variance Table

Model 1: logMaxAbund ~ logMass
Model 2: logMaxAbund ~ logMass + Diet
Res.Df RSS Df Sum of Sq F Pr(>F)
1 52 25.184
2 48 21.837 4 3.3477 1.8396 0.1366
'data.frame': 646 obs. of 15 variables:
$ abund : num 5 0.2 0.4 0 0 0 0 0 0 0 ...
$ Present : Factor w/ 2 levels "Absent","Present": 1 1 1 2 2 2 2 2 2 2 ...
$ clDD : num 5543 5750 5395 5920 6152 ...
$ clFD : num 83.5 67.5 79.5 66.7 57.6 59.2 59.5 51.5 47.4 46.3 ...
$ clTmi : num 9 9.6 8.6 11.9 11.6 10.8 10.8 11.6 13.6 13.5 ...
$ clTma : num 32.1 31.4 30.9 31.9 32.4 32.1 32.3 33 33.5 33.4 ...
$ clTmn : num 15.2 15.7 14.8 16.2 16.8 ...
$ clP : num 140 147 148 143 141 ...
$ NDVI : int -56 -44 -36 -49 -42 -49 -48 -50 -64 -58 ...
$ broadleaf: num 0.3866 0.9516 0.9905 0.0506 0.2296 ...
$ conif : num 0.0128 0.0484 0 0.9146 0.7013 ...
$ grass : num 0 0 0 0 0 0 0 0 0 0 ...
$ crop : num 0.2716 0 0 0.0285 0.044 ...
$ urban : num 0.2396 0 0 0 0.0157 ...
$ wetland : num 0 0 0 0 0 0 0 0 0 0 ...
Régression linéaire multiple dans R
Exécuter la régression multiple de l'abondance en fonction des variables clTma + NDVI + grass
summary(lm.mult)
Régression pas à pas
Méthode des moindres carrés
- la plus utilisée et correspond à la méthode par défaut dans R
valeur Yi observé (mesurée) à Xi
valeur Yi prédite à Xi
Moyenne de tout les Yi
VE:
résidus
VR:
variance régression
VT:
variance totale
R2
= VR/VT
Révision: Régression linéaire simple
Shapiro-Wilk normality test

data: resid(anov1)
W = 0.98, p-value = 0.4982
Bartlett test of homogeneity of variances

data: logMaxAbund by Diet
Bartlett's K-squared = 7.4728, df = 4, p-value = 0.1129
Exécution du modèle linéaire dans R
Exécution du modèle linéaire dans R
Hétéroscédastique
Non-linéaire
La fonction
abline()
permet d'ajouter la ligne
La fonction
hist()
permet de faire un histogramme
Suppositions non respectées - quelle est la cause?
Hypothèse
Essayez-les et comparez les sorties!
Rappel : ANCOVA
Solutions possibles

Exécuter une ANOVA dans R
Comparer les sorties
anov1 <- lm(logMaxAbund ~ Diet, data=bird)
anova(anov1)
aov1 <- aov(logMaxAbund ~ Diet, data=bird)
summary(aov1)
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
R2-aj: le modèle explique 11.28% de la variabilité de l'abondance de dickcissels
* Le partitionnement de la variation sera expliqué plus largement à l'Atelier 10.
Full transcript