Matériel requis

Cet atelier requiert les dernières versions de RStudio et de R.

Vous devez également avoir installé ces librairies R:

• Hmisc
• labdsv
• MASS
• vegan
• ggplot2

install.packages(c('Hmisc', 'labdsv', 'MASS', 'vegan', 'ggplot2'))

Vous pouvez télécharger les jeux de données sur r.qcbs.ca/fr/workshops/r-workshop-10.

Objectifs d'apprentissage

Objectif: Réaliser des analyses multivariées avancées sur des données de communauté.

En utilisant des méthodes d'ordination sous contrainte, nous apprendrons à:

• Explorer comment les variables environnementales expliquent des patrons de composition en espéces à travers différents sites.

• Décrire et prédire ces tendances entre la composition des communautés et des variables environnementales.

Introduction

Cet atelier est une continuation de l'atelier précédent (#9), qui a offert un aperçu des analyses multivariées non-contraintes:

• Mesures de distances et transformations de données
• Groupement hiérarchique
• Ordinations non-contraintes (PCA, PCoA, CA, nmDS)
  

Celles-ci permettent de relever des patrons dans la structure des communautés d'espèces ou des descripteurs.

L'atelier #10 montre des analyses permettant d'explorer comment les variables environnementales expliquent ces patrons.

Introduction

Cet atelier se concentrera sur les analyses sous contraintes:

• Analyse canonique de redondances (ACR ou RDA)
• Analyse canonique de redondances partielle (ACR partielle ou RDA partielle)
• Partitionnement de la variation
• Analyse linéaire discriminante (ALD ou LDA)
  
  

Ces analyses nous permettrons de décrire et de prédire les relations entre la structure des communautés et les variables environnementales.

On pourra alors tester des hypothèses!

Code et données

Lien vers le matériel de l'atelier: github.com/QCBSRworkshops/workshop10.

Téléchargez le code R et les données requises pour cet atelier:

• Code R
• Données:
  • DoubsEnv
  • DoubsSpe
  • DoubsSpa
  • Données test pour l'analyse linéaire discriminante

Pour télécharger les données, faire clic droit + Sauvegarder sur la page qui ouvre.

Librairies R

Assurez-vous de télécharger, d'installer et de charger les librairies suivants:

• Hmisc
• labdsv
• MASS
• vegan

Suivre l'atelier

Quelques conseils:

• Créez votre propre code (ou commentez le code R fourni);
• Évitez de copier-coller, ou d'exécuter le code directement du script fourni;
• N'oubliez pas de bien définir votre répertoire de travail vers le dossier contenant les fichiers requises pour l'atelier.

Introduction aux données

Rivière Doubs (Verneaux 1973)

Charger les données

Assurez-vous que les fichiers se trouvent dans votre répertoire de travail!

Chargez la matrice d'abondances d'espèces (doubsspe.csv).

# Assurez vous que les fichiers se trouvent dans votre répertoire de travail!
spe <- read.csv("data/doubsspe.csv", row.names = 1)
spe <- spe[-8,] # Supprimer site 8 (pas d'espèces).

Chargez la matrice de données environnementales (doubsenv.csv).

env <- read.csv("data/doubsenv.csv", row.names = 1)
env <- env[-8,] # Supprimer site 8 (pas d'espèces).

Note: N'exécuter qu'une seule fois!

Données d'abondances d'espèces

Explorons la matrice des abondances de poissons.

names(spe) # noms d'objets (espèces)
#  [1] "CHA" "TRU" "VAI" "LOC" "OMB" "BLA" "HOT" "TOX" "VAN" "CHE" "BAR" "SPI"
# [13] "GOU" "BRO" "PER" "BOU" "PSO" "ROT" "CAR" "TAN" "BCO" "PCH" "GRE" "GAR"
# [25] "BBO" "ABL" "ANG"
dim(spe) # dimensions de la matrice
# [1] 29 27

Et, si on veut plus de détails sur les espèces:

head(spe) # 5 premières lignes
str(spe) # structure d'objets de la matrice
summary(spe) # statistiques descriptives des objets (min, moyenne, max, etc.)

Distribution des abondances d'espèces

Explorons la structure de la communauté.

# Compter la fréquence d'espèces dans chaque classe d'abondance
ab <- table(unlist(spe))
# Visualiser cette distribution
barplot(ab, las = 1,
        xlab = "Abundance class", ylab = "Frequency",
        col = grey(5:0/5))

Note: Il y a beaucoup de zéros.

Distribution des abondances d'espèces

Combien y a-t-il de zéros?

sum(spe == 0)
# [1] 408

Quelle proportion de l'ensemble des données cela représente-t-il ?

sum(spe == 0)/(nrow(spe)*ncol(spe))
# [1] 0.5210728

Transformation des données d'abondances

Plus de 50% des données d'abondances sont en fait des absences. C'est élevé, mais pas inhabituel pour ce type de données.

Par contre, on ne veut pas que les double zéros amplifient artificiellement la similarité entre sites.

• Pour éviter ceci, on peut transformer les données d'abondances.

Nous appliquerons alors une transformation Hellinger:

# la fonction decostand() dans la libraire vegan nous facilite la tâche:
library(vegan)
spe.hel <- decostand(spe, method = "hellinger")

Données environnementales

Explorons les données environnementales.

names(env) # noms des objets (variables environnementales)
#  [1] "das" "alt" "pen" "deb" "pH"  "dur" "pho" "nit" "amm" "oxy" "dbo"
dim(env) # dimensions de la matrice
# [1] 29 11
head(env) # 5 premières lignes
#    das alt  pen  deb  pH dur  pho  nit  amm  oxy dbo
# 1  0.3 934 48.0 0.84 7.9  45 0.01 0.20 0.00 12.2 2.7
# 2  2.2 932  3.0 1.00 8.0  40 0.02 0.20 0.10 10.3 1.9
# 3 10.2 914  3.7 1.80 8.3  52 0.05 0.22 0.05 10.5 3.5
# 4 18.5 854  3.2 2.53 8.0  72 0.10 0.21 0.00 11.0 1.3
# 5 21.5 849  2.3 2.64 8.1  84 0.38 0.52 0.20  8.0 6.2
# 6 32.4 846  3.2 2.86 7.9  60 0.20 0.15 0.00 10.2 5.3

Et, si on veut plus de détails sur les variables environnementales:

str(env) # structure des objets
summary(env) # statistiques descriptives (min, moyenne, max, etc.)

Colinéarité

# On peut également détecter (visuellement) les colinéarités entres variables:
heatmap(abs(cor(env)), # corrélation de Pearson (note: ce sont des valeurs absolues!)
        col = rev(heat.colors(6)),
        Colv = NA, Rowv = NA)
legend("topright",
       title = "R de Pearson",
       legend =  round(seq(0,1, length.out = 6),1),
       y.intersp = 0.7, bty = "n",
       fill = rev(heat.colors(6)))

Standardisation des données

Il est impossible de comparer les effets de variables qui ont des unités différentes.

Avant d'effectuer les analyses qui suivent, les données doivent donc être standardisées.

# standardiser les données
env.z <- decostand(env, method = "standardize")
# centrer les données (moyenne ~ 0)
round(apply(env.z, 2, mean), 1)
# das alt pen deb  pH dur pho nit amm oxy dbo 
#   0   0   0   0   0   0   0   0   0   0   0
# réduire les données (écart type = 1)
apply(env.z, 2, sd)
# das alt pen deb  pH dur pho nit amm oxy dbo 
#   1   1   1   1   1   1   1   1   1   1   1

Analyses canoniques

Les analyses canoniques nous permettent:

• d'identifier les relations entre un ensemble de variables réponses et un ensemble de variables explicatives
• de tester des hypothèses écologiques à propos de ces relations
• de faire des prédictions

Analyse canonique de redondances (RDA)

L'analyse canonique de redondances est une ordination sous contraintes.

• La RDA est une extension directe de la régression multiple.
• Elle modélise l’effet d’une matrice explicative sur une matrice réponse (au lieu d'une seule variable réponse).

Les variables peuvent être quantitatives, qualitatives, ou binaires (0/1).

• transformez et standardisez les variables avant d'effectuer une RDA.

Effectuer une RDA dans R

Étape 1: Transformer et/ou standardiser les données.

# On utilisera nos données explicatives standardisées
# Enlever la variable "distance from the source" (colinéarité avec autres variables)
env.z <- subset(env.z, select = -das)

Étape 2: Effectuer une RDA.

# Modèlise l'effect de tous les variables environnementales sur la composition en espèces des communautés
spe.rda <- rda(spe.hel ~ ., data = env.z)

Étape 3: Extraire les résultats de la RDA.

summary(spe.rda)
...
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025     1.0000
# Constrained    0.3689     0.7341
# Unconstrained  0.1336     0.2659
...

RDA output in R

Pour interpréter les résultats d'une RDA, on peut d'abord se concentrer sur cette partie clé de la sortie:

...
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025     1.0000
# Constrained    0.3689     0.7341
# Unconstrained  0.1336     0.2659
...

• Constrained Proportion: variance de $Y$ expliquée par $X$ (73.41%)
• Unconstained Proportion: variance in $Y$ non expliquée par (26.59%)

Comment présenteriez-vous ces résultats?

• Les variables environnementales mesurées expliquent 73.41% de la variation dans la composition en espèces des communautés de poissons dans la rivière Doubs.

Sélection de variables

Une sélection progressive peut être effectuée afin de sélectionner les variables explicatives significatives.

Quelles variables contribuent de façon significative au pouvoir explicatif du modèle?

# Sélection progressive de variables:
fwd.sel <- ordiR2step(rda(spe.hel ~ 1, data = env.z), # modèle le plus simple
               scope = formula(spe.rda), # modèle "complet"
               direction = "forward",
               R2scope = TRUE, # limité par le R2 du modèle "complet"
               pstep = 1000,
               trace = FALSE) # mettre TRUE pour voir le processus du sélection!

Essentiellement, on ajoute une variable à la fois au modèle, et on retient la variable si elle augmente significativement le (R^2) ajusté du modèle.

Sélection de variables

• Quelles variables ont été sélectionnées?

fwd.sel$call
# rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z)

• Quel est le R2 ajusté d'une RDA incluant seulement les variables significatives?

# Écrire notre nouveau modèle
spe.rda.signif <- rda(spe.hel ~ alt + oxy + dbo, data = env.z)
# vérifier son R2 ajusté
RsquareAdj(spe.rda.signif)
# $r.squared
# [1] 0.5894243
# 
# $adj.r.squared
# [1] 0.5401552

Tester la significativité d'une RDA

Utilisez anova.cca() pour tester la significativité globale de notre RDA.

anova.cca(spe.rda.signif, permutations = 1000)
...
#          Df Variance      F   Pr(>F)    
# Model     3  0.29619 11.963 0.000999 ***
# Residual 25  0.20632                    
# ---
...

On peut aussi tester la significativité de chaque variable avec by = 'term'!

anova.cca(spe.rda.signif, permutations = 1000, by = "term")
...
# Model: rda(formula = spe.hel ~ alt + oxy + dbo, data = env.z)
#          Df Variance       F   Pr(>F)    
# alt       1 0.164856 19.9759 0.000999 ***
# oxy       1 0.082426  9.9877 0.000999 ***
# dbo       1 0.048909  5.9264 0.001998 ** 
# Residual 25 0.206319                     
...

Représentation graphique des RDAs

Une RDA permet la visualisation simultanée des variables réponses et explicatives (i.e. espèces et variables environnementales).

Comme pour la PCA, on doit choisir entre deux options de cadrage:

Triplot RDA: Cadrage de type 1

ordiplot(spe.rda.signif,
         scaling = 1,
         type = "text")

Triplot RDA: Cadrage de type 2

ordiplot(spe.rda.signif,
         scaling = 2,
         type = "text")

Configuration des triplots RDA

Les fonctions plot() et ordiplot() produisent des triplots rapidement et facilement, mais on peut aussi configurer les graphiques avec l'extraction de scores avec scores() et leur visualisation avec points(), text(), et arrows().

Voir le script et le livre de l'atelier 10 pour coder la figure!

Défi 1

Effectuez une RDA pour modéliser les effets des variables environnementales sur l’abondance des espèces d'acariens.

# Charger les données d'abondance des espèces d'acariens
data("mite")
# Charger les données environnementales
data("mite.env")

decostand()
rda()
ordiR2step()
anova.cca()
ordiplot()

Défi 1: Solution

Étape 1: Transformer et standardiser les données.

# Transformer les données d'abondances
mite.spe.hel <- decostand(mite, method = "hellinger")
# Standardiser les données environmentales quantiatives
mite.env$SubsDens <- decostand(mite.env$SubsDens, method = "standardize")
mite.env$WatrCont <- decostand(mite.env$WatrCont, method = "standardize")

Défi 1: Solution

Étape 2: Sélectionner les variables environnementales.

# RDA avec tous les variables environnementales
mite.spe.rda <- rda(mite.spe.hel ~ ., data = mite.env)
# Sélection progressive des variables environnementales significatives
fwd.sel <- ordiR2step(rda(mite.spe.hel ~ 1, data = mite.env),
                      scope = formula(mite.spe.rda),
                      direction = "forward",
                      R2scope = TRUE, pstep = 1000, trace = FALSE)
fwd.sel$call
# rda(formula = mite.spe.hel ~ WatrCont + Shrub + Substrate + Topo, 
#     data = mite.env)

Défi 1: Solution

Étape 3: Effectuer la RDA et extraire le R2 ajusté.

# Refaire la RDA avec seulement les variables significatives
mite.spe.rda.signif <- rda(mite.spe.hel ~ WatrCont + Shrub +
                           Substrate + Topo + SubsDens,
                           data = mite.env)
# Calculer le R2 ajusté
RsquareAdj(mite.spe.rda.signif)$adj.r.squared
# [1] 0.4367038

Défi 1: Solution

Étape 4: Tester la significativité globale du modèle.

anova.cca(mite.spe.rda.signif, step = 1000)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens, data = mite.env)
#          Df Variance     F Pr(>F)    
# Model    11  0.20759 5.863  0.001 ***
# Residual 58  0.18669                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Les variables environnementales sélectionnées expliquent 43.7% (p = 0.001) de la variation dans la composition de communautés des acariens.

Défi 1: Solution

Étape 5: Visualiser le modèle à l'aide d'un triplot!

ordiplot(mite.spe.rda.signif,
         scaling = 1,
         frame = F,
         main = "Cadrage 1")

ordiplot(mite.spe.rda.signif,
         scaling = 2,
         frame = F,
         main = "Cadrage 2")

RDA partielle

La RDA partielle est un cas particulier de la RDA qui permet de tenir compte de covariables.

En d'autres mots, on peut modéliser les effets linéaires de la matrice $X$ sur la matrice $Y$, tout en contrôlant pour l'effet de la matrice de covariables $W$.

Applications de la RDA partielle

Puisque la RDA partielle contrôle pour l'effet de covariables, on peut:

• Évaluer l’effet de variables environnementales sur la composition des communautés en prenant en compte des covariables de moindre intérêt.
• Isoler les effets d'un ou plusieurs groupes de variables explicatives.

Par exemple, la RDA partielle est souvent utilisée pour séparer les effets de variables environnementales des effets de variables spatiales.

Exemple: RDA partielle sur les données Doubs

Dans R, on peut faire une RDA partielle avec la fonction rda().

Par exemple, évaluons l'effet de la chimie de l'eau sur l’abondance des poissons (spe.hel) en tenant compte de covariables topographiques.

# Divisez le tableau de données environnementales en deux:
# variables topographiques et chimiques
env.topo <- subset(env.z, select = c(alt, pen, deb))
env.chem <- subset(env.z, select = c(pH, dur, pho, nit, amm, oxy, dbo))
# Faire la RDA partielle
spe.partial.rda <- rda(spe.hel, env.chem, env.topo)

spe.partial.rda <- rda(spe.hel ~ pH + dur + pho + nit + amm + oxy + dbo +
                       Condition(alt + pen + deb), # covariables ici
                       data = env.z)

Interprétation de la sortie d'une RDA partielle

summary(spe.partial.rda)

...
# Partitioning of variance:
#               Inertia Proportion
# Total          0.5025     1.0000
# Conditioned    0.2087     0.4153
# Constrained    0.1602     0.3189
# Unconstrained  0.1336     0.2659
...

• Conditioned Proportion: variance de $Y$ expliquée par $W$ (41.53%)
• Constrained Proportion: variance de $Y$ expliquée par $X$ (31.89%)
• Unconstained Proportion: variance de $Y$ non expliquée (26.59%)
  
  Comment présenteriez-vous ces résultats?

La chimie de l’eau explique 31.89% de l’abondance des espèces de poissons, tandis que la topographie explique 41.53% de la variation en abondances des poissons.

Interprétation d'une RDA partielle

Il nous manque encore quelques détails importants pour notre interprétation du modèle!

# Extraire le R2 ajusté du modèle
RsquareAdj(spe.partial.rda)$adj.r.squared
# [1] 0.2413464

# Évaluer la significativité statistique du modèle
anova.cca(spe.partial.rda, step = 1000)
...
# Permutation test for rda under reduced model
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem, Z = env.topo)
#          Df Variance      F Pr(>F)    
# Model     7  0.16024 3.0842  0.001 ***
# Residual 18  0.13360                  
...

Représentation graphique

On peut visualiser les effets des variables environnementales sur la communauté de poissons avec la fonction ordiplot().

ordiplot(spe.partial.rda,
         scaling = 2,
         main = "Rivière Doubs - Cadrage 2")

Note: Les variables topographiques ne sont pas représentées. Pourquoi?

Défi 2

rda()
summary()
RsquareAdj()
anova.cca() # voir l'argument 'by' dans ?anova.cca

Défi 2: Solution

Étape 1: Transformer et standardiser les données.

Nos données sont déjà transformés et standardisés!

Étape 2: Faire la RDA partielle:

mite.spe.subs <- rda(mite.spe.hel ~ Shrub + Topo
                     + Condition(SubsDens + WatrCont + Substrate),
                     data = mite.env)
# Extraire les résultats
summary(mite.spe.subs)
...
# Partitioning of variance:
#               Inertia Proportion
# Total         0.39428    1.00000
# Conditioned   0.16891    0.42839
# Constrained   0.03868    0.09811
# Unconstrained 0.18669    0.47350
...

Shrub et Topo expliquent 9.8% de la variation de l’abondance de mites, tandis que le substrat explique 42.8% de cette variation.

Défi 2: Solution

Étape 3: Interpréter les résultats!

• Quel pourcentage de la variance est expliqué par les variables environnementales?

RsquareAdj(mite.spe.subs)$adj.r.squared
# [1] 0.08327533

• Le modèle est-il significatif?

anova.cca(mite.spe.subs, step = 1000)
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ Shrub + Topo + Condition(SubsDens + WatrCont + Substrate), data = mite.env)
#          Df Variance     F Pr(>F)    
# Model     3 0.038683 4.006  0.001 ***
# Residual 58 0.186688                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Défi 2: Solution

• Quels axes sont significatifs?

anova.cca(mite.spe.subs, step = 1000, by = "axis")
# Permutation test for rda under reduced model
# Forward tests for axes
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(formula = mite.spe.hel ~ Shrub + Topo + Condition(SubsDens + WatrCont + Substrate), data = mite.env)
#          Df Variance      F Pr(>F)    
# RDA1      1 0.027236 8.4618  0.001 ***
# RDA2      1 0.008254 2.5643  0.022 *  
# RDA3      1 0.003193 0.9919  0.389    
# Residual 58 0.186688                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Partitionnement de la variation

Divise la variation d’une matrice de variables réponses en 2, 3, ou 4 matrices de variables explicatives.

• e.g. variables locales vs. à large échelle
• e.g. abiotique vs. biotique

Partitionnement de la variation

Partitionnement de la variation dans R

Pour démontrer le partitionnement de la variation dans R, nous allons partitionner la variation de la composition des espèces de poissons entre les variables chimiques et topographiques.

spe.part.all <- varpart(spe.hel, env.chem, env.topo)
spe.part.all$part # extraire résultats
# No. of explanatory tables: 2 
# Total variation (SS): 14.07 
#             Variance: 0.50251 
# No. of observations: 29 
# 
# Partition table:
#                      Df R.squared Adj.R.squared Testable
# [a+c] = X1            7   0.60579       0.47439     TRUE
# [b+c] = X2            3   0.41526       0.34509     TRUE
# [a+b+c] = X1+X2      10   0.73414       0.58644     TRUE
# Individual fractions                                    
# [a] = X1|X2           7                 0.24135     TRUE
# [b] = X2|X1           3                 0.11205     TRUE
# [c]                   0                 0.23304    FALSE
# [d] = Residuals                         0.41356    FALSE
# ---
# Use function 'rda' to test significance of fractions of interest

Diagramme de Venn

plot(spe.part.all,
     Xnames = c("Chem", "Topo"), # noms des matrices explicatives
     bg = c("seagreen3", "mediumpurple"), alpha = 80,
     digits = 2,
     cex = 1.5)

Tester la significativité

• La significativité de la fraction partagée [b] ne peut pas être testée.
• Mais, on peut tester la significativité des autres fractions!

Significativité: X1 [a+b]

anova.cca(rda(spe.hel, env.chem))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem)
#          Df Variance      F Pr(>F)    
# Model     7  0.30442 4.6102  0.001 ***
# Residual 21  0.19809                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significativité: X2 [b+c]

anova.cca(rda(spe.hel, env.topo))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.topo)
#          Df Variance     F Pr(>F)    
# Model     3  0.20867 5.918  0.001 ***
# Residual 25  0.29384                 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significativité: Fractions individuelles

anova.cca(rda(spe.hel, env.chem, env.topo))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.chem, Z = env.topo)
#          Df Variance      F Pr(>F)    
# Model     7  0.16024 3.0842  0.001 ***
# Residual 18  0.13360                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significativité: Fractions individuelles

anova.cca(rda(spe.hel, env.topo, env.chem))
# Permutation test for rda under reduced model
# Permutation: free
# Number of permutations: 999
# 
# Model: rda(X = spe.hel, Y = env.topo, Z = env.chem)
#          Df Variance      F Pr(>F)    
# Model     3 0.064495 2.8965  0.001 ***
# Residual 18 0.133599                  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Défi 3

data("mite.pcnm")

ordiR2step()
varpart()
anova.cca(rda())
plot()

Défi 3: Solution

Étape 1: Sélection de variables spatiales significatives

# Modèle RDA avec tous les variables spatiales
full.spat <- rda(mite.spe.hel ~ ., data = mite.pcnm)
# Sélection progressive des variables spatiales
spat.sel <- ordiR2step(rda(mite.spe.hel ~ 1, data = mite.pcnm),
               scope = formula(full.spat),
               R2scope = RsquareAdj(full.spat)$adj.r.squared,
               direction = "forward",
               trace = FALSE)
spat.sel$call
# rda(formula = mite.spe.hel ~ V2 + V3 + V8 + V1 + V6 + V4 + V9 + 
#     V16 + V7 + V20, data = mite.pcnm)

Défi 3: Solution

Étape 2: Créer des sous-groupes de variables explicatives.

# Variables de substrat
mite.subs <- subset(mite.env, select = c(SubsDens, WatrCont))
# Variables spatiales significatives
mite.spat <- subset(mite.pcnm,
                    select = names(spat.sel$terminfo$ordered))
                    # pour rapidement accèder aux variables sélectionnées

Défi 3: Solution

Étape 3: Partitionnement de la variation de la matrice d'abondances.

mite.part <- varpart(mite.spe.hel, mite.subs, mite.spat)
mite.part$part$indfract # extraire résultats
#                 Df R.squared Adj.R.squared Testable
# [a] = X1|X2      2        NA    0.05901929     TRUE
# [b] = X2|X1     10        NA    0.19415929     TRUE
# [c]              0        NA    0.24765221    FALSE
# [d] = Residuals NA        NA    0.49916921    FALSE

• Quelle est la proportion de variance expliquée par le substrat?

  • 5.9%

• Quelle est la proportion de variance expliquée par l'espace?

  • 19.4%

Défi 3: Solution

Étape 4: Quelles sont les fractions significatives?

[a]: Substrat seulement

anova.cca(rda(mite.spe.hel, mite.subs, mite.spat))
...
# Model: rda(X = mite.spe.hel, Y = mite.subs, Z = mite.spat)
#          Df Variance      F Pr(>F)    
# Model     2 0.025602 4.4879  0.001 ***
# Residual 57 0.162583                  
...

[c]: Espace seulement

anova.cca(rda(mite.spe.hel, mite.spat, mite.subs))
...
# Model: rda(X = mite.spe.hel, Y = mite.spat, Z = mite.subs)
#          Df Variance      F Pr(>F)    
# Model    10  0.10286 3.6061  0.001 ***
# Residual 57  0.16258                  
...

Défi 3: Solution

Étape 5: Visualiser les résultats avec un diagramme de Venn.

plot(mite.part,
     digits = 2,
     Xnames = c("Subs", "Space"), # titre des fractions
     cex = 1.5,
     bg = c("seagreen3", "mediumpurple"), # ajoutez des couleurs!
     alpha = 80)

Défi 3: Solution

Alors, quels sont les effets du substrat et de l'espace sur les abondances d'espèces de mites?

En groupes, résumez et interprétez ces résultats comme si vous écriviez un article. Indice: Pourquoi trouve-t-on un effet si important de l'espace?

Analyse discriminante linéaire (LDA)

Une LDA divise une variable réponse en groupes selon un facteur en trouvant une combinaison de variables qui donne la meilleure séparation possible entre les groupes. Ceci est utile parce que:

• La LDA détermine si une matrice de variables indépendantes explique bien un groupement établi a priori;

• Et permet donc de prédire la classification de nouvelles données!

  e.g. prédire l’appartenance d’une espèce de poisson à un groupe (marin vs. eau douce) selon sa morphologie

LDA dans R: Rivière Doubs

Généralement, les variables environnementales changent avec la latitude. On peut donc poser la question suivante:

LDA dans R: Rivière Doubs

Commençons par charger les données spatiales des sites:

# charger les données spatiales des sites Doubs:
spa <- read.csv("data/doubsspa.csv", row.names = 1)
spa$site <- 1:nrow(spa) # assigner un chiffre par site
spa <- spa[-8,] # enlever le site #8

Ensuite, on peut classifier les sites dans 3 groupes de latitudes:

spa$group <- NA # créer colonne "group"
spa$group[which(spa$y < 82)] <- 1
spa$group[which(spa$y > 82 & spa$y < 156)] <- 2
spa$group[which(spa$y > 156)] <- 3

LDA dans R: Rivière Doubs

Visualisons ces regroupements par latitude:

LDA dans R

Note: Normalement, nous devons vérifier que les matrices de covariance des variables explicatives sont homogènes (voir Borcard et al. 2011). Au besoin, on peut utiliser la fonction betadisper() du paquet vegan.

Pour les besoins de l'atelier, on passera directement à la LDA:

# charger la libraire requise
library(MASS)
# faire la LDA
LDA <- lda(env, spa$group)

LDA dans R

Nos sites sont maintenant réorganisés en groupes qui sont les plus distincts possibles, à partir des variables environnementales.

# prédire les groupes à partir de la LDA
lda.plotdf <- data.frame(group = spa$group, lda = predict(LDA)$x)

LDA dans R: Exactitude du groupement

On peut ensuite voir comment les sites sont classifiés en groupes, et si cette classification est exacte.

# classification des objets en fonction de la LDA
spe.class <- predict(LDA)$class
# probabilités que les objets appartiennent à chaque groupe a posteriori
spe.post <- predict(LDA)$posterior
# tableau des classifications a priori et prédites
(spe.table <- table(spa$group, spe.class))
#    spe.class
#      1  2  3
#   1  7  0  0
#   2  0  6  0
#   3  0  0 16
# proportion de classification correcte
diag(prop.table(spe.table, 1))
# 1 2 3 
# 1 1 1

Tous les sites ont été correctement classifiés dans leur groupe de latitude en fonction des variables environnementales.

LDA dans R: Prédictions

On peut maintenant utiliser la LDA pour classifier de nouveaux sites dans les groupes de latitude.

Essayons de prédire la classification de cinq nouveaux sites à l'aide de la LDA:

# charger les nouvelles données
classify.me <- read.csv("data/classifyme.csv", header = TRUE)
# enlever das
classify.me <- subset(classify.me, select = -das)
# prédire le groupement des nouvelles données
predict.group <- predict(LDA, newdata = classify.me)
# prédire la classification pour chaque site
predict.group$class
# [1] 1 1 1 3 3
# Levels: 1 2 3

Défi 4

data(mite.xy)

lda()
predict()
table()
diag()

Défi 4: Solution

Étape 1: Créer quatre groupes de latitude avec des étendues égales.

# numéroter les sites
mite.xy$site <- 1:nrow(mite.xy)
# trouver une étendue égale de latitudes par groupe
(max(mite.xy[,2])-min(mite.xy[,2]))/4
# [1] 2.4
# classifier les sites dans 4 groupes de latitude
mite.xy$group <- NA # nouvelle colonne "group"
mite.xy$group[which(mite.xy$y < 2.5)] <- 1
mite.xy$group[which(mite.xy$y >= 2.5 & mite.xy$y < 4.9)] <- 2
mite.xy$group[which(mite.xy$y >= 4.9 & mite.xy$y < 7.3)] <- 3
mite.xy$group[which(mite.xy$y >= 7.3)] <- 4

Défi 4: Solution

Étape 2: Faire la LDA.

LDA.mite <- lda(mite.env[,1:2], mite.xy$group)

Étape 3: Vérifier l'exactitude du groupement.

# classification des objects en fonction de la LDA
mite.class <- predict(LDA.mite)$class
# tableeau de classifications  (prior versus predicted)
(mite.table <- table(mite.xy$group, mite.class))
#    mite.class
#      1  2  3  4
#   1  9  4  2  0
#   2  2 11  4  0
#   3  1  2 14  2
#   4  0  0  3 16
# proportion de classifications exactes
diag(prop.table(mite.table, 1))
#         1         2         3         4 
# 0.6000000 0.6470588 0.7368421 0.8421053

Défi 4: Solution

On peut maintenant répondre à la question du défi avec cette partie du code:

# proportion of correct classification
diag(prop.table(mite.table, 1))
#         1         2         3         4 
# 0.6000000 0.6470588 0.7368421 0.8421053

Alors, quelle proportion de sites ont été classifiés correctement dans le groupe 1? Le groupe 2?

• 60% des sites ont été classifiés correctement dans le groupe 1, et 64.7% dans le groupe 2.

C'est ça! On a réussi!