class: center, middle, inverse, title-slide # Atelier 4: Modèles linéaires ## Série d’ateliers R du CSBQ ### Centre de la Science de la Biodiversité du Québec --- --- class: inverse, center, middle # Concepts importants ## Définir la moyenne et la variation --- ## Moyenne La moyenne est une mesure de la valeur moyenne d'une population (*x*): `$$\bar{x} = \frac{1}{N} \sum_{i=1}^{n} x_{i}$$` --- ## Variation - La variation est la dispersion des observations autour de la moyenne - Écart moyen - Variance - Écart type - Coefficient de variation **Mais qu'est-ce que c'est l'écart ?** `$$D_{i} = |x_{i} - \bar{x}|$$` --- ## Variation **L'écart**: `$$D_{i} = |x_{i} - \bar{x}|$$` -- **L'écart moyen**: `$$D = \frac{1}{N} \sum_{i=1}^{n} |x_{i} - \bar{x}|$$` -- Au lieu de valeurs absolues, nous pouvons également mettre la valeur au carré, donnant la **variance**: `$$V = \frac{1}{N} \sum_{i=1}^{n} {(x_{i} - \bar{x})}^2$$` --- ## Variation Mais en mettant chaque valeur au carré, ces variables ne sont plus en unités significatives On fait donc la racine carrée de la **variance** ( `\(V\)` ), donnant l'**écart type**: `$$\sigma = \sqrt{V}$$` -- L'écart type relatif, en pourcentage, est le **coefficient de variation**: `$$cv = \frac{\sigma}{\bar{x}}$$` --- class: inverse, center, middle # Les modèles linéaires --- ## Les modèles linéaires Relation linéaire entre une variable réponse ( `\(Y\)` ) et une/des explicative ( `\(X\)` ), en utilisant les concepts de **moyenne** et **variance** - `\(Y\)` : variable que vous voulez expliquer (une seule variable réponse) - `\(X\)` : variable pour expliquer variable réponse (une ou plusieurs variables explicatives) - `\(Y\)` : doit être quantitative - `\(X\)` : quantitative ou qualitative - `\(\epsilon\)` : ce qui n'est pas expliqué par la ou les variables explicatives  résidus ou erreur --- ## Définir des modèles linéaires Rassembler tous ensemble : `$$Y_{i} = \beta_{0} + \beta_{1} x_{i1} + \cdots + \beta_{p} x_{ip} + \epsilon_{i}$$` - `\(Y_i\)` est la variable réponse - `\(β_0\)` est l'ordonnée à l'origine de la droite de régression - `\(β_1\)`est le coefficient de variation de la `\(1^{ère}\)` variable explicative - `\(β_p\)` est le coefficient de variation de la `\(p^{ème}\)` variable explicative - `\(x_{i1}\)` est la variable explicative quantitative pour la `\(1^{ère}\)` observation - `\(x_ip\)` est la variable explicative quantitative pour la `\(p^{ème}\)` observation - `\(ε_i\)` sont les résidus du modèle (i.e. la variance inexpliquée) --- ## Le but des modèles linéaires - Le but d'un modèle linéaire est de trouver la meilleure estimation des paramètres (les variables `\(\beta\)`), puis d'évaluer la qualité de l'ajustement du modèle - Plusieurs méthodes ont été développées pour calculer l'intercept et les coefficient de modèles linéaires - Le choix approprié dépend du type de variables explicatives considérées et de leur nombre .center[.large[Le concept général de ces méthodes consiste à minimiser les résidus]] --- ## Objectif d'enseignement .center[  ] --- ## Conditions de base du modèle linéaire 1. Les résidus sont indépendants 2. Les résidus suivent une distribution normale 3. Les résidus ont une moyenne de 0 4. Les résidus sont homoscédastiques (i.e. leur variance est constante) .alert[Ces 4 conditions concernent les résidus, et non les variables réponses ou explicatives] .small[.comment[Dans les sections suivantes, nous ne répétons pas les conditions ci-dessus pour chaque modèle parce que ces conditions de base s'appliquent à tous les modèles linéaires]] --- ## Flux de travail .center[  ] - Visualiser les données - Créer un modèle - Tester les 4 conditions de base du modèle - Ajuster le modèle si les conditions de base ne sont pas respectées - Interpréter les résultats du modèle --- class: inverse, center, middle # Régression linéaire simple --- ## Régression linéaire simple - Type de modèle linéaire qui contient seulement une variable explicative continue `$$Y_i = \beta_0 + \beta_1 x_i + \epsilon_i$$` - Estimation de l'**ordonnée à l'origine** ( `\(\beta_0\)` ) et un **coefficient de corrélation** ( `\(\beta_1\)` ) - Méthode des moindres carrés - méthode la plus couramment utilisée (défaut sur R) --- ## Méthode des moindres carrés .pull-left[ .center[] ] .pull-right[ **Suppositions** - `\(Y_i\)` : valeur observé (mesurée) pour `\(X_i\)` - `\(\widehat{Y}_i\)` : valeur prédite pour `\(X_i\)` - `\(\bar{Y}\)` : moyenne de tout les `\(Y_i\)` - `\(V_E\)` : résidus (erreur) - `\(V_R\)` : variance expliqué par la régression - `\(V_T\)` : variance totale - `\(R^2 = \frac{V_R}{V_T}\)` ] --- ## Effectuer un modèle linéaire .small[ **Étape 1**. Exécuter votre modèle linéaire **Étape 2**. Vérifier les suppositions ] .pull-left[.center[]] .pull-right[.center[]] .pull-left[.center[*Suppositions sont satisfaites ?*] .small[**Étape 3**. Estimer les paramètres de régression, test de signification, tracer votre modèle ]] .pull-right[.center[*Suppositions non satisfaites ?*] .small[*Pouvez-vous transformer vos variables (est-ce justifié) ?*] .pull-left[.center[]] .pull-right[.center[]] .small[ .pull-left[ Oui: retourner à l'étape 1 avec des variables transformées ] .pull-right[ Non: essayer GLM qui pourrait mieux convenir aux données ]]] --- ## Exécution du modèle linéaire dans R **Étape 1**. créer votre modèle linéaire Dans R, la fonction `lm()` est utilisée pour ajuster un modèle linéaire ```r lm1 <- lm(Y~X) ``` - `lm1` : Nouvel objet contenant le modèle linéaire - `Y` : Variable réponse - `X` : Variable indépendante --- ## Exécution du modèle linéaire dans R Télécharger les donées <span style="color:blue"> *birdsdiet* </span>: ```r bird <- read.csv("birdsdiet.csv") ``` Visualisez le tableau de la structure des données en utilisant la fonction `str()` : ```r str(bird) # 'data.frame': 54 obs. of 7 variables: # $ 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 ... ``` --- ## Exécution du modèle linéaire dans R Variable réponse : **abondance d'oiseaux**  num : quantitative Variable explicative : **masse**  num : quantitative ```r str(bird) # 'data.frame': 54 obs. of 7 variables: # $ 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 ... ``` Nous voulons d'abord vérifier si l'abondance maximale des oiseaux (`maxAbund`) est une fonction de la masse des oiseaux (`Mass`) ```r lm1 <- lm(MaxAbund ~ Mass, data = bird) ``` --- ## Exécution du modèle linéaire dans R **Étape 2**. Vérifier les suppositions avec les graphiques diagnostics ```r opar <- par(mfrow=c(2,2)) plot(lm1) ``` - `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 La sortie comprend les quatre graphiques diagnostics de la fonction `lm()` --- ## Graph. #1 - Résidus vs valeurs prédites Exemple d'indépendance .comment[(ce que nous recherchons !)] - Devrait montrer une dispersion de points sans patron <img src="workshop04-fr_files/figure-html/unnamed-chunk-10-1.png" width="396" style="display: block; margin: auto;" /> --- ## Graph. #1 - Résidus vs valeurs prédites Example de non-indépendance .comment[(ce que nous ne voulons pas !)] <img src="workshop04-fr_files/figure-html/unnamed-chunk-11-1.png" width="612" style="display: block; margin: auto;" /> - 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.) --- ## Graph. #2 - échelle localité Devrait montrer une dispersion de points sans patron <img src="workshop04-fr_files/figure-html/unnamed-chunk-12-1.png" width="576" style="display: block; margin: auto;" /> .pull-left[.center[]] .pull-right[.center[]] .pull-left[.center[Aucun patron dans les résidus]] .pull-right[.center[Forte tendance dans les résidus]] --- ## Graph. # 3 - Normal QQ - Compare la distribution (quantiles) des résidus aux quantiles d'une distribution normale - Si les points se situent de façon linéaire sur la ligne 1: 1, les résidus suivent une distribution normale <img src="workshop04-fr_files/figure-html/unnamed-chunk-13-1.png" width="576" style="display: block; margin: auto;" /> .pull-left[.center[]] .pull-right[.center[]] .pull-left[.center[C'est bien !]] .pull-right[.center[Pas très bien...]] --- ## Graph. # 4 - Résidus vs effet de levier - 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** - Les points de levier avec une forte influence peuvent être identifiés avec une **distance de Cook supérieure à 0,5** --- ## Effet levier vs influence <img src="workshop04-fr_files/figure-html/unnamed-chunk-14-1.png" width="360" style="display: block; margin: auto;" /> --- ## Effet levier vs influence <img src="workshop04-fr_files/figure-html/unnamed-chunk-15-1.png" width="648" style="display: block; margin: auto;" /> .pull-left[.center[]] .pull-right[.center[]] .pull-left[.center[Aucune valeur influente]] .pull-right[.center[Effet de levier élevé et influence raisonnable]] <br /> <br /> <br /> .comment[Ici, le point 29 a un effet de levier élevé, mais son influence est acceptable (à l'intérieur des limites de la distance Cook de 0,5)] --- ## Effet levier vs influence Effet de levier et influence élevée <br /> .pull-left[ <img src="workshop04-fr_files/figure-html/unnamed-chunk-16-1.png" width="288" style="display: block; margin: auto;" /> ] .pull-right[ - Points en dehors de la limite de 0,5 de la distance Cook - Ces points ont trop d'influence sur la régression ] .alert[Vous ne devriez jamais supprimer les valeurs aberrantes si vous n'avez pas de bonnes raisons de le faire (ex: erreur de mesure)] --- ## **Étape 2**. Vérifier les suppositions de `lm1` ```r par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) plot(lm1) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-17-1.png" width="396" style="display: block; margin: auto;" /> --- ## Suppositions non-respectées - quelle est la cause? Traçons le graphique Y ~ X avec la droite de régression et les histogrammes de Y et X pour explorer leurs distributions .small[ ```r plot(bird$MaxAbund ~ bird$Mass) abline(lm1) # adds the best-fit line hist(bird$MaxAbund) # hist() produces a histogram of the variable hist(bird$Mass) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-18-1.png" width="648" style="display: block; margin: auto;" /> ] --- ## Suppositions non-respectées - quelle est la cause? Vérifions la normalité des données à l'aide d'un test de `Shapiro-Wilk` et d'un test d'asymétrie (**`skewness`**) : ```r shapiro.test(bird$MaxAbund) # # Shapiro-Wilk normality test # # data: bird$MaxAbund # W = 0.5831, p-value = 3.872e-11 shapiro.test(bird$Mass) # # Shapiro-Wilk normality test # # data: bird$Mass # W = 0.54667, p-value = 1.155e-11 ``` .comment[Dans les deux cas, les distributions ne sont pas normales ] --- ## Suppositions non-respectées - quelle est la cause? Vérifions la normalité des données à l'aide d'un test de `Shapiro-Wilk` et d'un test d'asymétrie (**`skewness`**) : ```r skewness(bird$MaxAbund) # [1] 3.167159 skewness(bird$Mass) # [1] 3.146062 ``` .comment[La valeur positive indique que la distribution des données est décalée vers la gauche] --- ## Transformer les données - Normalisons les données en appliquant une transformation `log10()` - Ajoutons ces variables transformées à notre base de données ```r bird$logMaxAbund <- log10(bird$MaxAbund) bird$logMass <- log10(bird$Mass) ``` **Étape 1**: Exécuter une régression linéaire sur les données transformées ```r lm2 <- lm(logMaxAbund ~ logMass, data = bird) ``` --- ## **Étape 2**: Vérifier les suppositions de `lm2` ```r par(mfrow=c(2,2), mar=c(3,4,1.15,1.2)) plot(lm2) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-23-1.png" width="424.8" style="display: block; margin: auto;" /> .comment[.center[Beaucoup mieux !]] --- ## **Étape 2**: Vérifier les suppositions de `lm2` ```r plot(logMaxAbund ~ logMass, data=bird) abline(lm2) hist(log10(bird$MaxAbund)) hist(log10(bird$Mass)) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-24-1.png" width="648" style="display: block; margin: auto;" /> --- ## **Étape 3**: Estimer les paramètres et leur seuil de signification La fonction `summary()` est utilisée pour obtenir les paramètres, leur importance, etc .small[ ```r summary(lm2) # # Call: # lm(formula = logMaxAbund ~ logMass, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.93562 -0.39982 0.05487 0.40625 1.61469 # # 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 ``` ] --- ## **Étape 3**: Estimer les paramètres et leur seuil de signification Nous pouvons aussi extraire les paramètres du modèle, par exemple : ```r lm2$coef # (Intercept) logMass # 1.6723673 -0.2361498 summary(lm2)$coefficients # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.6723673 0.2471519 6.766557 1.166186e-08 # logMass -0.2361498 0.1169836 -2.018658 4.869342e-02 summary(lm2)$r.squared # [1] 0.07267022 ``` --- ## Discussion de groupe - 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` ? .small[ ```r summary(lm2) # # Call: # lm(formula = logMaxAbund ~ logMass, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.93562 -0.39982 0.05487 0.40625 1.61469 # # 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 ``` ] --- ## Discussion de groupe Pouvons-nous améliorer le modèle si nous analysons que les oiseaux terrestres? .comment[Vous pouvez exclure des objets en utilisant `=!`] ```r lm3 <- lm(logMaxAbund~logMass, data=bird, subset=!bird$Aquatic) # removes aquatic birds (= TRUE) # or equivalently lm3 <- lm(logMaxAbund~logMass, data=bird, subset=bird$Aquatic == 0) ``` ```r # Examine the diagnostic plots par(mfrow=c(2,2)) plot(lm3) summary(lm3) # Compare both models par(mfrow=c(1,2)) plot(logMaxAbund~logMass, data=bird) plot(logMaxAbund~logMass, data=bird, subset=!bird$Aquatic) ``` --- ## Plot `R2-adj` a changé de 0,05 à 0,25 quand nous avons exclu les oiseaux aquatiques : ```r par(mfrow=c(1,2), mar = c(4, 4, 3, 1)) plot(logMaxAbund~logMass, data=bird, main = 'Tous les oiseaux') abline(lm2, col = 'red') plot(logMaxAbund~logMass, data=bird, subset=!bird$Aquatic, main = 'Oiseaux terrestres') abline(lm3, col = 'red') ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-30-1.png" width="504" style="display: block; margin: auto;" /> --- ## Défi 1 ![:cube]() - Examiner la relation entre `log(MaxAbund)` et `log(Mass)` chez les passereaux ("passerine birds") - Sauvegarder l'objet du modèle sous `lm4` .comment[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] - Comparer la variance expliquée par `lm2`, `lm3` and `lm4` --- ## Défi 1 - Solution ![:cube]() <br> ```r # Run the model lm4 <- lm(logMaxAbund ~ logMass, data=bird, subset=bird$Passerine == 1) summary(lm4) # # Call: # lm(formula = logMaxAbund ~ logMass, data = bird, subset = bird$Passerine == # 1) # # Residuals: # Min 1Q Median 3Q Max # -1.24644 -0.20937 0.02494 0.25192 0.93624 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.2429 0.4163 2.985 0.00661 ** # logMass 0.2107 0.3076 0.685 0.50010 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # Residual standard error: 0.5151 on 23 degrees of freedom # Multiple R-squared: 0.02, Adjusted R-squared: -0.02261 # F-statistic: 0.4694 on 1 and 23 DF, p-value: 0.5001 ``` --- ## Défi 1 - Solution ![:cube]() ```r # diagnostic plots par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-32-1.png" width="468" style="display: block; margin: auto;" /> --- ## Défi 1 - Solution ![:cube]() Comparer la variance expliquée par `lm2`, `lm3` and `lm4` ```r # Recall: we want adj.r.squared summary(lm2)$adj.r.squared # [1] 0.05483696 summary(lm3)$adj.r.squared # [1] 0.2484544 summary(lm4)$adj.r.squared # [1] -0.02260731 ``` .comment[Le meilleur modèle parmi les trois est `lm3` *(seulement les oiseaux terrestres)*] --- ## Objectif d'enseignement .center[  ] --- class: inverse, center, middle # ANOVA ## Test-t ## ANOVA à un critère de classification ## ANOVA à deux critères de classification --- ## ANOVA Variable réponse continue **Variables explicatives catégoriques** - Deux niveaux ou plus (groupes) .large[.center[Compare la variation intra-groupe et inter-groupe afin de déterminer si les moyennes des groupes diffèrent]] Somme des carrés : variance intra-traitement *vs* variance inter-traitement Si variance inter traitements `\(>\)` variance intra traitements: - la variable explicative a un effet plus important que l'erreur aléatoire - variable explicative est donc susceptible d'influencer significativement la variable réponse --- ## Types d'ANOVA 1. ANOVA à un critère de classification - Une variable explicative catégorique avec au moins 2 niveaux - S'il y a 2 niveaux, un **test t** peut être utilisé alternativement 2. ANOVA à deux critères de classification - Deux variables explicatives catégoriques ou plus - Chaque facteur peut avoir plusieurs niveaux - Les interactions entre chaque variable explicative catégorique doivent être testées 3. Mesures répétées ? - L'ANOVA peut être utilisée pour des mesures répétées, mais ce sujet n'est pas abordé dans cet atelier - Modèle linéaire mixte peut également être utilisé pour ce type de données (voir l'atelier 6) --- class: inverse, center, middle # Test T --- ## Test T - **Variable réponse**  quantitative - **Variable explicative**  qualitative avec **2 niveaux** **Suppositions** - Les résidus suivent une distribution normale - Les variances des groupes sont homogènes .comment[Le test est plus robuste lorsque la taille de l'échantillon est plus élevée et lorsque les groupes ont des tailles égales] --- ## Exécuter un test T dans R Vous pouvez utiliser la fonction `t.test()` ```r t.test(Y~X2, data= data, alternative = "two.sided") ``` - `Y`: variable réponse - `X2`: facteur (2 niveaux) - `data`: nom du jeu de doné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 Vous pouvez donc aussi utiliser la fonction `lm()` ```r lm.t <-lm(Y~X2, data = data) anova(lm.t) ``` --- ## Exécuter un test T dans R .large[Les oiseaux aquatiques sont-ils plus lourds que les oiseaux terrestres ?] - Variable réponse : `Bird mass`  num: continue - Variable explicative : `Aquatic`  2 niveaux : 1 ou 0 (oui ou non) --- ## Exécuter un test T dans R Premièrement, visualiser les données à l'aide de la fonction `boxplot()` ```r boxplot(logMass ~ Aquatic, data = bird, names = c("Non aquatique", "Aquatique")) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-34-1.png" width="468" style="display: block; margin: auto;" /> --- ## Exécuter un test T dans R Testons l'homogénéité des variances avec la fonction `var.test()` ```r var.test(logMass ~ Aquatic, data = bird) # # 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 ``` .comment[Le rapport des variances n'est pas statistiquement différent de 1, celles-ci peuvent donc être considérées comme égales] .comment[Nous pouvons maintenant procéder au test t !] --- ## Exécuter un test T dans R ```r ttest1 <- t.test(logMass ~ Aquatic, var.equal = TRUE, data = bird) # Or use lm() ttest.lm1 <- lm(logMass ~ Aquatic, data=bird) ``` .comment[Spécifie que l'homogénéité des variances est respectée] Vérifiez que `t.test()` et `lm()` donnent le même modèle : ```r ttest1$statistic^2 # t # 60.3845 anova(ttest.lm1)$F # [1] 60.3845 NA # réponse : F=60.3845 dans les deux cas ``` .comment[Lorsque la supposition d'égalité de variance est confirmée, t^2 = F] --- ## Exécuter un test T dans R 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 ```r 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 in group 0 mean in group 1 # 1.583437 2.908289 ``` .small[.comment[Il existe une différence entre la masse des oiseaux aquatiques et terrestres - `p-value`]] .small[.comment[Regardez les moyennes des 2 groupes]] --- ## Non respect des suppositions - **Correction de Welch** : lorsque les écarts entre les groupes ne sont pas égaux (par défaut dans R !) - **Test de Mann-Whitney** : l'équivalent **non paramétrique** du test de t lorsque les suppositions ne sont pas respectées - **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) --- ## Discussion de groupe .large[Les oiseaux aquatiques sont-ils plus lourds que les oiseaux terrestres ?] ```r # Unilateral t-test uni.ttest1 <- t.test(logMass ~ Aquatic, var.equal = TRUE, data = bird, alternative = "less") ``` .comment[Qu'avez-vous conclu ?] --- ## Discussion de groupe ```r uni.ttest1 # # 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 ``` .small[ Oui, les oiseaux aquatiques sont plus lourds que les oiseaux terrestres : p-value = 0.0000000001468087 ] --- class: inverse, center, middle # ANOVA --- ## Analyse de Variance (ANOVA) Généralisation du test t à `\(>2\)` groupes, et/ou ≥ `\(2\)` facteurs explicatifs Décomposition de la variance observée de la variable réponse en effets additifs d'un ou de plusieurs facteurs et de leurs interactions <br> `$$Y = \underbrace{\mu}_{\Large{\text{moyenne globale de la variable réponse}\atop\text{sur tous les individus}}} + \overbrace{\tau_{i}}^{\Large{\text{Le résultat moyen sur}\atop\text{tous les individus du groupe i}}} + \underbrace{\epsilon}_{\text{Résidus}}$$` --- ## 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** --- ## Exécuter une ANOVA dans R ##### 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 ```r str(bird) # 'data.frame': 54 obs. of 9 variables: # $ 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 ... # $ logMaxAbund: num 0.475 1.577 2.383 0.643 0.656 ... # $ logMass : num 2.855 0.724 1.554 2.077 2.499 ... ``` --- ## Visualiser les données Visualisons tout d'abord les données avec la fonction `boxplot()` ```r boxplot(logMaxAbund ~ Diet, data = bird, ylab = expression("log"[10]*"(Abondance maximale)"), xlab = 'Régime alimentaire') ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-42-1.png" width="504" style="display: block; margin: auto;" /> --- ## Visualiser les données 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()` ```r med <- sort(tapply(bird$logMaxAbund, bird$Diet, median)) boxplot(logMaxAbund ~ factor(Diet, levels = names(med)), data = bird, ylab = expression("log"[10]*"(Abondance maximale)"), xlab = 'Régime alimentaire') ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-43-1.png" width="504" style="display: block; margin: auto;" /> --- ## Visualiser les données Une autre façon de visualiser graphiquement les tailles d’effet est d’utiliser la fonction `plot.design()` .small[ ```r plot.design(logMaxAbund ~ Diet, data=bird, ylab = expression("log"[10]*"(Abondance maximale)")) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-44-1.png" width="432" style="display: block; margin: auto;" /> ] .comment[Les niveaux d'un facteur le long d'une ligne verticale, et la valeur globale de la réponse dans une ligne horizontale] --- ## Exécuter une ANOVA à un critère de classification dans R Il est de nouveau possible d'utiliser la fonction `lm()` ```r anov1 <- lm(logMaxAbund ~ Diet, data = bird) ``` Il y a aussi une fontion spécifique pour l'analyse de la variance dans R `aov()` ```r aov1 <- aov(logMaxAbund ~ Diet, data = bird) ``` .comment[Essayez-les et comparez les sorties !] --- ## Exécuter une ANOVA à un critère de classification dans R Comparer les sorties .small[ ```r anova(anov1) # Analysis of Variance Table # # Response: logMaxAbund # Df Sum Sq Mean Sq F value Pr(>F) # Diet 4 5.1059 1.27647 2.8363 0.0341 * # Residuals 49 22.0521 0.45004 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r 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 ``` ] --- ## Vérifier les suppositions **Test de Bartlett**: égalité de la variance entre les groupes .small[ ```r bartlett.test(logMaxAbund ~ Diet, data = bird) # # Bartlett test of homogeneity of variances # # data: logMaxAbund by Diet # Bartlett's K-squared = 7.4728, df = 4, p-value = 0.1129 ``` ] **Test de Shapiro-Wilk**: normalité des résidus .small[ ```r shapiro.test(resid(anov1)) # # Shapiro-Wilk normality test # # data: resid(anov1) # W = 0.97995, p-value = 0.4982 ``` ] .comment[Les deux tests sont non-significatifs; les résidus du modèle peuvent être considérés normaux et les variances homogènes] --- ## 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 ```r data$logY <- log10(data$Y) ``` * Voir le wiki de l'atelier 1 pour les règles de transformation de données * 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*) ```r kruskal.test(Y~X, data) ``` --- ## Sorties de notre modèle ANOVA Triage en ordre alphabétique des niveaux et comparaison au niveau de référence (`Insect`) .small[ ```r summary(anov1) # # Call: # lm(formula = logMaxAbund ~ Diet, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.85286 -0.32972 -0.08808 0.47375 1.56075 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.1539 0.1500 7.692 5.66e-10 *** # DietInsectVert -0.6434 0.4975 -1.293 0.2020 # DietPlant 0.2410 0.4975 0.484 0.6303 # DietPlantInsect 0.4223 0.2180 1.938 0.0585 . # DietVertebrate -0.3070 0.2450 -1.253 0.2161 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # Residual standard error: 0.6709 on 49 degrees of freedom # Multiple R-squared: 0.188, Adjusted R-squared: 0.1217 # F-statistic: 2.836 on 4 and 49 DF, p-value: 0.0341 ``` ] --- ## Sorties de notre modèle ANOVA D'autre part, si nous utilisons `lm()` .pull-left2[ .small[ ```r summary.lm(aov1) # # Call: # aov(formula = logMaxAbund ~ Diet, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.85286 -0.32972 -0.08808 0.47375 1.56075 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.1539 0.1500 7.692 5.66e-10 *** # DietInsectVert -0.6434 0.4975 -1.293 0.2020 # DietPlant 0.2410 0.4975 0.484 0.6303 # DietPlantInsect 0.4223 0.2180 1.938 0.0585 . # DietVertebrate -0.3070 0.2450 -1.253 0.2161 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # Residual standard error: 0.6709 on 49 degrees of freedom # Multiple R-squared: 0.188, Adjusted R-squared: 0.1217 # F-statistic: 2.836 on 4 and 49 DF, p-value: 0.0341 ``` ]] .pull-right2[ <br><br> .comment[Différence significative entre les groupes, mais nous ne savons pas lesquels !]] --- ## Test a posteriori Lorsque l'ANOVA détecte un effet significatif de la variable explicative, un test post-hoc avec la fonction `TkeyHSD()`, doit être effectué pour déterminer quel(s) tratement(s) diffère(nt) .pull-left2[ .small[ ```r TukeyHSD(aov(anov1), ordered = TRUE) # 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 ``` ]] .pull-right2[ <br><br> .comment[Seuls `Vertebrate` et `PlantInsect` diffèrent] ] --- ## Représentation graphique Représentation graphique de l'ANOVA à l'aide de la fonction `barplot()` .small[ ```r sd <- tapply(bird$logMaxAbund, bird$Diet, sd) means <- tapply(bird$logMaxAbund, bird$Diet, mean) n <- length(bird$logMaxAbund) se <- 1.96*sd/sqrt(n) bp <- barplot(means, ylim = c(0, max(bird$logMaxAbund) - 0.5)) epsilon = 0.1 segments(bp, means - se, bp, means + se, lwd=2) # barres verticales segments(bp - epsilon, means - se, bp + epsilon, means - se, lwd = 2) # barres horizontales segments(bp - epsilon, means + se, bp + epsilon, means + se, lwd = 2) # barres horizontales ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-56-1.png" width="504" style="display: block; margin: auto;" /> ] --- class: inverse, center, middle # ANOVA à deux critères de classification --- ## ANOVA à deux critères de classification Plus d'un facteur - ANOVA avec un facteur: `aov <- lm(Y ~ X, data)` - ANOVA avec deux ou plus facteurs: `aov <- lm(Y ~ X * Z * ..., data)` .comment[lorsque vous utilisez le symbole "*" avec `lm()`, le modèle inclut les effets de chaque facteur séparément, ainsi que leur interaction] .comment[lorsque vous utilisez le symbole "+" avec `lm()`, le modèle inclut les effets de chaque facteur séparément (pas d'interaction)] `aov <- lm(Y ~ X + Z + ..., data)` --- ## ANOVA à deux critères de classification Exemple d'interaction non significative .small[ ```r aov <- lm(Y ~ X * Z, data) summary(aov) # 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 ``` 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 ```r aov <- lm(Y ~ X + Z, data) ``` ] --- ## Défi 2 ![:cube]() Testez si l'abondance maximale `log(MaxAbund)` varie à la fois en fonction du régime alimentaire (`Diet`) et de l'habitat (`Aquatic`) - .comment[INDICE: Examinez les facteurs Diet, Aquatic et leur interaction avec une ANOVA à deux critères de classificatioN e.g. `lm(Y ~ A*B)`] - .comment[où A est le premier facteur, B le deuxième et "*" décrit l'interaction] --- ## Défi 2 - Solution ![:cube]() .tiny[ ```r anov2 <- lm(logMaxAbund ~ Diet*Aquatic, data = bird) summary(anov2) # # Call: # lm(formula = logMaxAbund ~ Diet * Aquatic, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.9508 -0.2447 0.0000 0.3584 1.1558 # # Coefficients: (1 not defined because of singularities) # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.23447 0.17324 7.126 6.64e-09 *** # DietInsectVert -0.86989 0.67097 -1.296 0.2014 # DietPlant 0.16043 0.49001 0.327 0.7449 # DietPlantInsect 0.35358 0.23395 1.511 0.1377 # DietVertebrate -0.95449 0.33772 -2.826 0.0070 ** # Aquatic -0.26858 0.31630 -0.849 0.4003 # DietInsectVert:Aquatic 0.56034 0.96976 0.578 0.5663 # DietPlant:Aquatic NA NA NA NA # DietPlantInsect:Aquatic 0.05516 0.73821 0.075 0.9408 # DietVertebrate:Aquatic 1.24044 0.49408 2.511 0.0157 * # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # Residual standard error: 0.6482 on 45 degrees of freedom # Multiple R-squared: 0.3038, Adjusted R-squared: 0.18 # F-statistic: 2.454 on 8 and 45 DF, p-value: 0.02688 ``` ] --- ## Défi 2 - Solution ![:cube]() .small[ ```r anov2 <- lm(logMaxAbund ~ Diet*Aquatic, data = bird) 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 ``` ] .comment[Le seul terme significatif du modèle est le facteur régime alimentaire] --- ## Objectif d'enseignement .center[  ] --- class: inverse, center, middle # ANCOVA --- ## Analyse de covariance (ANCOVA) - Combinaison de l'ANOVA et de la régression linéaire - Les variables explicatives sont un mélange de variables quantitatives (covariable) et qualitatives (facteurs) `$$Y = \mu + \text{Effets principaux des facteurs} + \\ \text{Interactions entre facteurs} + \\ \text{Effets principaux des covariables} + \\ \text{Interactions entre covariables et facteurs} + \epsilon$$` --- ## Rappel : ANCOVA En plus des suppositions des modèles linéaires, les modèles **ANCOVA** doivent respecter : - Les covariables ont toutes la **même étendue de valeurs** - Les variables sont **fixes** - Les variables catégoriques et continues sont **indépendantes** <br> .small[ .comment[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)]] .small[.comment[*Voir l'atelier 6 sur les modèles linéaires mixtes*]] --- ## Types d'ANCOVA 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 <br> ANCOVA fréquemment utilisées 1. **Une covariable et un facteur** 2. Une covariable et deux facteurs 3. Deux covariables et un facteur .small[.comment[Nous ne considérerons que le premier cas aujourd'hui, mais les deux autres sont similaires]] --- ## ANCOVA avec 1 covariable et 1 facteur Objectifs de l'analyse : 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 <br> .center[.alert[Si vous avez une interaction significative entre votre facteur et votre covariable, vous ne pouvez pas atteindre ces objectifs !]] --- ## ANCOVA avec 1 covariable et 1 facteur <br> <img src="workshop04-fr_files/figure-html/unnamed-chunk-59-1.png" width="720" style="display: block; margin: auto;" /> .pull-left2[.center[.pull-left[] .pull-right[]]] .pull-right2[.center[]] .center[.pull-left2[.small[Si l'interaction est significative, vous aurez un scénario qui ressemble à ceci]]] .pull-right2[.small[Si votre covariable et votre facteur sont significatifs, vous avez un cas comme celui-ci]] --- ## Comparez ANCOVA - moyennes ajustées 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 .small[ ```r ancova.exemple <- lm(Y ~ X*Z, data=data) # X = quantitative; Z = qualitative library(effects) adj.means.ex <- effect('Z', ancova.exemple) plot(adj.means.ex) ``` ] <img src="workshop04-fr_files/figure-html/unnamed-chunk-61-1.png" width="360" style="display: block; margin: auto;" /> --- ## ANCOVA avec 1 covariable et 1 facteur - Si seulement votre facteur est significatif, éliminer la covariable -> vous avez une **ANOVA** - Si seulement votre covariable est significative, éliminer le facteur -> vous avez une **régression linéaire simple** - Si votre interaction covariable * facteur est significative, vous voudrez peut-être tester quel(s) niveau(x) du facteur a(ont) des pentes différentes .alert[Vérifier vos suppositions ! ] - Très similaire à ce que vous avez fait précédemment --- ## Exécuter une ANCOVA dans R ##### L'abondance maximale varie-t-elle en fonction du régime alimentaire et la masse des oiseaux ? Variable réponse : **MaxAbund**  num : quantitative continue Variable explicatives : - **Diet**  facteur à 5 niveaux - **Mass**  numérique continue .small[ ```r str(bird) # 'data.frame': 54 obs. of 9 variables: # $ 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 ... # $ logMaxAbund: num 0.475 1.577 2.383 0.643 0.656 ... # $ logMass : num 2.855 0.724 1.554 2.077 2.499 ... ``` ] --- ## Défi 3 ![:cube]() 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`) ```r ancova.exemple <- lm(Y~X*Z, data=data) summary(ancova.exemple) ``` 2- Vérifiez si votre interaction est significative ```r ancova.exemple2 <- lm(Y~X+Z, data=data) summary(ancova.exemple2) ``` --- ## Défi 3 - Solution ![:cube]() <br> ```r 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 ``` .comment[Interaction entre `logMass` et `Diet` est non-significative] --- ## Défi 3 - Solution ![:cube]() Éliminer le terme d'interaction, puis ré-évaluer le modèle contenant les effets simples de `logMass` et `Diet` ```r ancov2 <- lm(logMaxAbund ~ logMass + Diet, data = bird) anova(ancov2) # Analysis of Variance Table # # Response: logMaxAbund # Df Sum Sq Mean Sq F value Pr(>F) # logMass 1 1.9736 1.97357 4.3382 0.04262 * # Diet 4 3.3477 0.83691 1.8396 0.13664 # Residuals 48 21.8367 0.45493 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## Objectif d'enseignement .center[  ] --- class: inverse, center, middle # Régression linéaire multiple --- ## Régression linéaire multiple - **Variables explicatives**  2 ou plusieurs variables continues - **Variable réponse**  1 variable continue Il s'agit d'une généralisation de la régression linéaire simple `$$Y_i = \alpha + \beta_1x_{1i}+\beta_2x_{2i}+\beta_3x_{3i}+...+\beta_kx_{ki} + \epsilon$$` Suppositions (.small[*En plus des suppositions habituelles des modèles linéaires*]) - La relation entre chaque variable explicative et la variable réponse est de type linéaire - Les variables explicatives sont orthogonales (pas de colinéarité) Si colinéarité 1. Gardez seulement **une** des variables colinéaires 2. Essayez une analyse multidimensionelle (voir l'atelier 9) 3. Essayez une analyse pseudo-orthogonale --- ## Régression linéaire multiple dans R 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`) .small[ ```r Dickcissel = read.csv("data/dickcissel.csv") str(Dickcissel) # '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 ... ``` ] --- ## Vérifier les suppositions La colinéarité : - Vérifier la colinéarité de toutes les variables explicatives et d'intérêt .small[ .pull-left[ ```r # select variables var <- c('clTma', 'NDVI', 'grass', 'abund') plot(Dickcissel[, var]) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-68-1.png" width="504" style="display: block; margin: auto;" /> ]] .pull-right[ <br> .small[.comment[ 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 ]]] --- ## Régression linéaire multiple dans R Exécuter la régression multiple de l'abondance (`abund`) en fonction des variables `clTma + NDVI + grass` ```r lm.mult <- lm(abund ~ clTma + NDVI + grass, data = Dickcissel) summary(lm.mult) ``` Vérifiez les autres suppositions, comme pour la régression linéaire simple ```r par(mfrow = c(2, 2)) plot(lm.mult) ``` --- ## Régression linéaire multiple dans R Exécuter la régression multiple de l'abondance (`abund`) en fonction des variables `clTma + NDVI + grass` .small[ ```r lm.mult <- lm(abund ~ clTma + NDVI + grass, data = Dickcissel) summary(lm.mult) # # Call: # 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 ``` ] --- ## Régression linéaire multiple dans R Vérifiez les autres suppositions, comme pour la régression linéaire simple ```r par(mfrow = c(2, 2)) plot(lm.mult) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-72-1.png" width="432" style="display: block; margin: auto;" /> --- ## Quel est le meilleur modèle ? .small[ Souvenez-vous du principe de parcimonie: expliquer le plus de variation avec le plus petit nombre de termes dans votre modèle  enlevez la variable qui est la moins significative ] .pull-left2[ .small[ ```r summary(lm.mult) # # Call: # 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 ``` ]] .pull-right2[ .small[.comment[ <br> Les 3 variables sont importantes. On garde tout ! R2-aj: le modèle explique 11.28% de la variabilité de l'abondance de dickcissels ]]] --- ## Quel est le meilleur modèle ? Il est important de noter que la variable réponse ne varie pas de façon linéaire avec les variables explicatives ```r plot(abund ~ clTma, data = Dickcissel) plot(abund ~ NDVI, data = Dickcissel) plot(abund ~ grass, data = Dickcissel) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-74-1.png" width="648" style="display: block; margin: auto;" /> .comment[Voir la **section avancée** sur la **régression polynomiale** pour la solution !] --- class: inverse, center, middle # Optionnel ## *si le temps le permet* --- ## Optionnel 1. Régression pas à pas 2. Interprétation des contrastes 3. ANOVA non équilibrée 4. Régression polynomiale 5. Partitionnement de la variation --- class: inverse, center, middle # Régression pas à pas --- ## Régression pas à pas Exécuter un modèle avec tout dedans sauf les variables de "Présent/absence" 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 .small[ ```r lm.full <- lm(abund ~ . - Present, data = Dickcissel) lm.step <- step(lm.full) ``` ] --- ## Régression pas à pas .pull-left[ .tiny[ ```r summary(lm.full) # # Call: # lm(formula = abund ~ . - Present, data = Dickcissel) # # Residuals: # Min 1Q Median 3Q Max # -32.151 -9.847 -2.864 4.215 170.774 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) -444.39158 35.48446 -12.524 < 2e-16 *** # clDD 0.35577 0.33576 1.060 0.28973 # clFD 1.06739 0.17534 6.088 1.99e-09 *** # clTmi -6.73898 0.76207 -8.843 < 2e-16 *** # clTma 3.40939 1.34342 2.538 0.01139 * # clTmn -110.10517 122.81542 -0.897 0.37032 # clP 0.14624 0.05080 2.879 0.00413 ** # NDVI 0.01949 0.05446 0.358 0.72048 # broadleaf 1.38425 3.67225 0.377 0.70634 # conif 0.65936 4.02652 0.164 0.86998 # grass 10.66584 5.18745 2.056 0.04018 * # crop -0.86060 3.35285 -0.257 0.79751 # urban -31.01974 22.69841 -1.367 0.17224 # wetland 34.06486 806.29941 0.042 0.96631 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # # Residual standard error: 19.9 on 632 degrees of freedom # Multiple R-squared: 0.3248, Adjusted R-squared: 0.3109 # F-statistic: 23.39 on 13 and 632 DF, p-value: < 2.2e-16 ``` ]] .pull-right[ Variables sélectionées par `step()` .tiny[ ```r summary(lm.step) # # 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 ``` ]] .small[.comment[Le modèle explique maintenant 31,44% de la variabilité de l'abondance de Dickcissel]] --- class: inverse, center, middle # Interprétation des contrastes --- ## Interprétation des contrastes .small[ Par défaut, `contr.treatment` compare chaque niveau du facteur à un niveau de référence L'estimation de l'ordonnée à l'origine 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 .comment[*Que remarquez-vous ?*] ```r tapply(bird$logMaxAbund, bird$Diet, mean) # Insect InsectVert Plant PlantInsect Vertebrate # 1.1538937 0.5104603 1.3948941 1.5761940 0.8468898 coef(anov1) # (Intercept) DietInsectVert DietPlant DietPlantInsect # 1.1538937 -0.6434334 0.2410004 0.4223003 # DietVertebrate # -0.3070039 coef(anov1)[1] + coef(anov1)[2] # InsectVert # (Intercept) # 0.5104603 coef(anov1)[1] + coef(anov1)[3] # Plant # (Intercept) # 1.394894 ``` ] --- ## Interprétation des contrastes Il se peut que vous vouliez définir un niveau de référence différent 1. Comparez le niveau `Plant` à tous les autres niveaux du facteur `Diet` ```r bird$Diet2 <- relevel(bird$Diet, ref="Plant") anov2 <- lm(logMaxAbund ~ Diet2, data = bird) summary(anov2) anova(anov2) ``` 2. Ordonner les niveaux selon leur médiane ```r bird$Diet2 <- factor(bird$Diet, levels=names(med)) anov2 <- lm(logMaxAbund ~ Diet2, data = bird) summary(anov2) anova(anov2) ``` .comment[Observez-vous un changement quant aux niveaux du facteur `Diet` qui sont significatifs ?] --- ## Interprétation des contrastes .comment[Un point important à remarquer à propos du contraste par défaut dans R (`contr.treatment`) est qu'il n'est PAS orthogonal] Pour être orthogonal : - Pour être orthogonal, les propriétés suivantes doivent être respectées: - La somme du produit de deux colonnes égale 0 ```r sum(contrasts(bird$Diet)[,1]) # [1] 1 sum(contrasts(bird$Diet)[,1]*contrasts(bird$Diet)[,2]) # [1] 0 ``` --- ## Interprétation des contrastes Changez les contrastes pour mettre les niveaux orthogonaux .small[ ```r options(contrasts=c("contr.helmert", "contr.poly")) sum(contrasts(bird$Diet)[,1]) # [1] 0 sum(contrasts(bird$Diet)[,1]*contrasts(bird$Diet)[,2]) # [1] 0 ``` ] .pull-left[ .tiny[ ```r anov3 <- lm(logMaxAbund ~ Diet, data = bird) summary(anov3) # # Call: # lm(formula = logMaxAbund ~ Diet, data = bird) # # Residuals: # Min 1Q Median 3Q Max # -1.85286 -0.32972 -0.08808 0.47375 1.56075 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 1.09647 0.14629 7.495 1.14e-09 *** # Diet1 -0.32172 0.24876 -1.293 0.2020 # Diet2 0.18757 0.17854 1.051 0.2986 # Diet3 0.13911 0.06960 1.999 0.0512 . # Diet4 -0.06239 0.05238 -1.191 0.2393 # --- # 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 ``` ]] .pull-right[ .small[Les contrastes Helmert vont contraster le deuxième niveau avec le premier, le troisième avec la moyenne des deux premiers niveaux, etc.] ] --- class: inverse, center, middle # ANOVA non équilibrée --- ## 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) ```r table(bird$Aquatic) # # 0 1 # 39 15 ``` L'ordre des covariables a changé la valeur de la somme des carrés ```r unb.anov1 <- lm(logMaxAbund ~ Aquatic + Diet, data = bird) unb.anov2 <- lm(logMaxAbund ~ Diet + Aquatic, data = bird) ``` --- ## ANOVA non équilibrée ```r anova(unb.anov1) # Analysis of Variance Table # # Response: logMaxAbund # Df Sum Sq Mean Sq F value Pr(>F) # Aquatic 1 0.2316 0.23157 0.5114 0.47798 # Diet 4 5.1926 1.29816 2.8671 0.03291 * # Residuals 48 21.7337 0.45279 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ```r anova(unb.anov2) # Analysis of Variance Table # # Response: logMaxAbund # Df Sum Sq Mean Sq F value Pr(>F) # Diet 4 5.1059 1.27647 2.8191 0.03517 * # Aquatic 1 0.3183 0.31834 0.7031 0.40591 # Residuals 48 21.7337 0.45279 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` --- ## ANOVA non équilibrée Maintenant essayez une `Anova()` de type III .pull-left[ .small[ ```r car::Anova(unb.anov1, type = "III") # Anova Table (Type III tests) # # Response: logMaxAbund # Sum Sq Df F value Pr(>F) # (Intercept) 18.9349 1 41.8186 4.837e-08 *** # Aquatic 0.3183 1 0.7031 0.40591 # Diet 5.1926 4 2.8671 0.03291 * # Residuals 21.7337 48 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ]] .pull-right[ .small[ ```r car::Anova(unb.anov2, type = "III") # Anova Table (Type III tests) # # Response: logMaxAbund # Sum Sq Df F value Pr(>F) # (Intercept) 18.9349 1 41.8186 4.837e-08 *** # Diet 5.1926 4 2.8671 0.03291 * # Aquatic 0.3183 1 0.7031 0.40591 # Residuals 21.7337 48 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ]] .comment[Que remarquez-vous en utilisant `Anova()` ?] --- class: inverse, center, middle # Régression polynomiale --- ## Régression polynomiale 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 - Un modèle polynômial ressemble à ceci : .center[$$\underbrace{2x^4}+\underbrace{3x}-\underbrace{2}$$] .comment[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 <br> .center[*Nous avons ici un polynôme de degré 4*] `$$2x^\overbrace{4} + 3x - 2$$` --- ## Régression polynomiale Lorsque vous connaissez le degré, vous pouvez lui donner un nom : <table> <thead> <tr> <th style="text-align:right;"> degre </th> <th style="text-align:left;"> Nom </th> <th style="text-align:left;"> Example </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;"> 0 </td> <td style="text-align:left;"> Constante </td> <td style="text-align:left;"> \(3\) </td> </tr> <tr> <td style="text-align:right;"> 1 </td> <td style="text-align:left;"> Linéaire </td> <td style="text-align:left;"> \(x+9\) </td> </tr> <tr> <td style="text-align:right;"> 2 </td> <td style="text-align:left;"> Quadratique </td> <td style="text-align:left;"> \(x^2-x+4\) </td> </tr> <tr> <td style="text-align:right;"> 3 </td> <td style="text-align:left;"> Cubique </td> <td style="text-align:left;"> \(x^3-x^2+5\) </td> </tr> <tr> <td style="text-align:right;"> 4 </td> <td style="text-align:left;"> Quartique </td> <td style="text-align:left;"> \(6x^4-x^3+x-2\) </td> </tr> <tr> <td style="text-align:right;"> 5 </td> <td style="text-align:left;"> Quintique </td> <td style="text-align:left;"> \(x^5-3x^3+x^2+8\) </td> </tr> </tbody> </table> --- ## Régression polynomiale 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) : ```r 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 polynomiale Comparez les modèles polynômiaux; .comment[quel modèle niché nous devrions sélectionner ?] Exécutez un résumé de ce modèle .tiny[ ```r print(summ_lm.linear) # [1] "Coefficients:" # [2] " Estimate Std. Error t value Pr(>|t|) " # [3] "(Intercept) 1.864566 2.757554 0.676 0.49918 " # [4] "clDD 0.001870 0.000588 3.180 0.00154 **" # [5] "Multiple R-squared: 0.01546,\tAdjusted R-squared: 0.01393 " # [6] "F-statistic: 10.11 on 1 and 644 DF, p-value: 0.001545" ``` ```r print(summ_lm.quad) # [1] "Coefficients:" # [2] " Estimate Std. Error t value Pr(>|t|) " # [3] "(Intercept) -1.968e+01 5.954e+00 -3.306 0.001 ** " # [4] "clDD 1.297e-02 2.788e-03 4.651 4.00e-06 ***" # [5] "I(clDD^2) -1.246e-06 3.061e-07 -4.070 5.28e-05 ***" # [6] "Multiple R-squared: 0.04018,\tAdjusted R-squared: 0.0372 " # [7] "F-statistic: 13.46 on 2 and 643 DF, p-value: 1.876e-06" ``` ```r print(summ_lm.cubic) # [1] "Coefficients:" # [2] " Estimate Std. Error t value Pr(>|t|)" # [3] "(Intercept) -1.465e+01 1.206e+01 -1.215 0.225" # [4] "clDD 8.612e-03 9.493e-03 0.907 0.365" # [5] "I(clDD^2) -1.628e-07 2.277e-06 -0.071 0.943" # [6] "I(clDD^3) -8.063e-11 1.680e-10 -0.480 0.631" # [7] "Multiple R-squared: 0.04053,\tAdjusted R-squared: 0.03605 " # [8] "F-statistic: 9.04 on 3 and 642 DF, p-value: 7.202e-06" ``` ] --- class: inverse, center, middle # Partitionnement de la variation --- ## Partitionnement de la variation Certaines variables explicatives de la **régression linéaire multiple** étaient fortement corrélées (c.-à-d. multicolinéarité) La colinéarité entre variables explicatives peut être détectée à l'aide de critères d'inflation de la variance (fonction `vif()` du packet `car`) - Les valeurs supérieures à 5 sont considérées colinéaires ```r mod <- lm(clDD ~ clFD + clTmi + clTma + clP + grass, data = Dickcissel) car::vif(mod) # clFD clTmi clTma clP grass # 13.605855 9.566169 4.811837 3.196599 1.165775 ``` --- ## Partitionnement de la variation .small[ 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)] .pull-left2[ .tiny[ ```r library(vegan) part.lm = varpart(Dickcissel$abund, Dickcissel[ ,c("clDD","clFD","clTmi","clTma","clP")], Dickcissel[ ,c("broadleaf","conif","grass","crop", "urban","wetland")]) part.lm # # Partition of variance in RDA # # Call: varpart(Y = Dickcissel$abund, X = Dickcissel[, c("clDD", # "clFD", "clTmi", "clTma", "clP")], Dickcissel[, c("broadleaf", # "conif", "grass", "crop", "urban", "wetland")]) # # Explanatory tables: # X1: Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")] # X2: Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")] # # No. of explanatory tables: 2 # Total variation (SS): 370770 # Variance: 574.84 # No. of observations: 646 # # Partition table: # Df R.squared Adj.R.squared Testable # [a+b] = X1 5 0.31414 0.30878 TRUE # [b+c] = X2 6 0.03654 0.02749 TRUE # [a+b+c] = X1+X2 11 0.32378 0.31205 TRUE # Individual fractions # [a] = X1|X2 5 0.28456 TRUE # [b] 0 0.02423 FALSE # [c] = X2|X1 6 0.00327 TRUE # [d] = Residuals 0.68795 FALSE # --- # Use function 'rda' to test significance of fractions of interest ``` ]] .pull-right2[ <br><br> .small[.comment[**Note** : les variables colinéaires n'ont pas besoin d'être enlevées avant l'analyse]] ] --- ## Partitionnement de la variation .pull-left[ .small[ ```r showvarparts(2) ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-98-1.png" width="432" style="display: block; margin: auto;" /> ```r ?showvarparts # With two explanatory tables, the fractions # explained uniquely by each of the two tables # are ‘[a]’ and ‘[c]’, and their joint effect # is ‘[b]’ following Borcard et al. (1992). ``` ]] .pull-right[ .small[ ```r plot(part.lm, digits = 2, bg = rgb(48,225,210,80, maxColorValue=225), col = "turquoise4") ``` <img src="workshop04-fr_files/figure-html/unnamed-chunk-100-1.png" width="432" style="display: block; margin: auto;" /> ]] .small[.comment[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%]] --- ## Partitionnement de la variation Tester si chaque fraction est significative - Climat seul ```r out.1 = rda(Dickcissel$abund, Dickcissel[ ,c("clDD", "clFD","clTmi","clTma","clP")], Dickcissel[ ,c("broadleaf","conif","grass","crop", "urban","wetland")]) ``` - Couverture du paysage seul ```r out.2 = rda(Dickcissel$abund, Dickcissel[ ,c("broadleaf","conif","grass","crop", "urban", "wetland")], Dickcissel[ ,c("clDD","clFD","clTmi", "clTma","clP")]) ``` --- ## Partitionnement de la variation .pull-left[ .small[ ```r # Climat seul anova(out.1, step = 1000, perm.max = 1000) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(X = Dickcissel$abund, Y = Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")], Z = Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")]) # Df Variance F Pr(>F) # Model 5 165.12 53.862 0.001 *** # Residual 634 388.72 # --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` ]] .pull-right[ .small[ ```r # Couverture du paysage seul anova(out.2, step = 1000, perm.max = 1000) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(X = Dickcissel$abund, Y = Dickcissel[, c("broadleaf", "conif", "grass", "crop", "urban", "wetland")], Z = Dickcissel[, c("clDD", "clFD", "clTmi", "clTma", "clP")]) # Df Variance F Pr(>F) # Model 6 5.54 1.5063 0.167 # Residual 634 388.72 ``` ]] .comment[ 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] --- class: inverse, center, bottom # Merci d'avoir participé !