9 Association, partitionnement et ordination

️ Objectifs spécifiques:

À la fin de ce chapitre, vous

  • serez en mesure d’effectuer des calculs permettant de mesurer des différence entre des observations, des groupes d’observation ou des variables observées
  • serez en mesure d’effection des analyses de partitionnement hiérarchiques et non-hiérarchiques
  • serez en mesure d’effectuer des calculs d’ordination à l’aide des techniques de réduction d’axe communes: analyse en composante principale, l’analyse de correspondance, l’analyse en coordonnées principales, analyse discriminante linéaire, l’analyse de redondance et l’analyse canonique des correspondances.

Les données écologiques incluent généralement plusieurs variables qui doivent être analysées conjointement. Les techniques pour l’analyse multivariée de données écologiques ont grandi en nombre et en complexité, laissant émerger l’écologie numérique, un nouveau domaine d’étude scientifique initié par Pierre Legendre et Louis Legendre dont l’ouvrage Numerical Ecology, aujourd’hui à sa troisième édition, reste un incontournable pour qui s’intéresse aux mathématiques sous-jacentes au domaine. Pour la rédaction de ces notes, c’est toutefois le livre Numerical ecology with R, écrit par Borcard et al. (2011) pour offrir un guide à qui voudrait une approche plus appliquée.

L’écologie numérique sera effleurée dans ce chapitre, qui introduit à trois concepts.

  1. Les associations permettent de quantifier la ressemblance ou la différence entre deux observation (échantillons) ou variables (descripteurs). Lorsque l’on a plus de deux variables ou plus de deux site, nous obtenons des matrices d’association.
  2. Le partitionnement permet de regrouper des observations ou des variables selon des métriques d’association.
  3. L’ordination vise par l’intermédiaire de techniques de réduction d’axe à mettre de l’ordre dans des données dont le nombre élevé de variables peut amener à des difficultés d’appréciation et d’interprétaion.

9.1 Espaces d’analyse

9.1.1 Abondance et occurence

L’abondance est le décompte d’espèces observées, tandis que l’occurence est la présence ou l’absence d’une espèce. Le tableau suivant contient des données d’abondance.

Ce tableau peut être rapidement transformé en données d’occurence, qui ne comprennent que l’information booléenne de présence (noté 1) et d’absence (noté 0).

L’espace des espèces (ou des variables ou descripteurs) est celui où les espèces forment les axes et où les sites sont positionnés dans cet espace. Il s’agit d’une perspective en mode R, qui permet principalement d’identifier quels espèces se retrouvent plus courrament ensemble.

## # A tibble: 4 x 3
##   `Bruant chanteur` Chardonneret `Mésange à tête noire`
##               <dbl>        <dbl>                  <dbl>
## 1                 1            0                     20
## 2                 0            9                      1
## 3                 5            6                      1
## 4                 2            0                      0

Dans l’espace des sites (ou les échantillons ou objets), on transpose la matrice d’abondance. On passe ici en mode Q, où chaque point est une espèce, et où l’on peut observer quels échantillons sont similaires.

##                           [,1] [,2] [,3] [,4]
## Bruant familier              1    0    0    3
## Citelle à poitrine rousse    1    0    0    0
## Colibri à gorge rubis        0    1    0    0
## Geai bleu                    3    2    0    0
## Bruant chanteur              1    0    5    2
## Chardonneret                 0    9    6    0
## Bruant à gorge blanche       1    0    0    0
## Mésange à tête noire        20    1    1    0
## Jaseur boréal               66    0    0    0

9.1.2 Environnement

L’espace de l’environnement comprend souvent un autre tableau contenant l’information sur l’environnement où se trouve les espèces: les coordonnées et l’élévation, la pente, le pH du sol, la pluviométrie, etc.

9.2 Analyse d’association

Nous utiliserons le terme association comme une mesure pour quantifier la ressemblance ou la différence entre deux objets (échantillons) ou variables (descripteurs).

Alors que la corrélation et la covariance sont des mesures d’association entre des variables (analyse en mode R), la similarité et la distance sont deux types de une mesure d’association entre des objets (analyse en mode Q). Une distance de 0 est mesurée chez deux objets identiques. La distance augmente au fur et à mesure que les objets sont dissociés. Une similarité ayant une valeur de 0 indique aucune association, tandis qu’une valeur de 1 indique une association parfaite. À l’opposé, la dissimilarité est égale à 1-similarité.

La distance peut être liée à la similarité par la relation:

\[distance=\sqrt{1-similarité}\]

ou

\[distance=\sqrt{dissimilarité}\]

La racine carrée permet, pour certains indices de similarité, d’obtenir des propriétés euclédiennes. Pour plus de détails, voyez le tableau 7.2 de Legendre et Legendre (2012).

Les matrices d’association sont généralement présentées comme des matrices carrées, dont les dimensions sont égales au nombre d’objets (mode Q) ou de vrariables (mode R) dans le tableau. Chaque élément (“cellule”) de la matrice est un indice d’association entre un objet (ou une variable) et un autre. Ainsi, la diagonale de la matrice est un vecteur nul (distance ou dissimilarité) ou unitaire (similarité), car elle correspond à l’association entre un objet et lui-même.

Puisque l’association entre A et B est la même qu’entre B et A, et puisque la diagonale retourne une valeur convenue, il est possible d’exprimer une matrice d’association en mode “compact”, sous forme de vecteur. Le vecteur d’association entre des objets A, B et C contiendra toute l’information nécessaire en un vecteur de trois chiffres, [AB, AC, BC], plutôt qu’une matrice de dimension \(3 \times 3\). L’impact sur la mémoire vive peut être considérable pour les calculs comprenant de nombreuses dimensions.

En R, les calculs de similarité et de distances peuvent être effectués avec le module vegan. La fonction vegdist permet de calculer les indices d’association en forme carrée.

Nous verons plus tard les méthodes de mesure de similarité et de distance plus loin. Pour l’instant, utilisons la méthode de Jaccard pour une démonstration sur des données d’occurence.

##           1         2         3         4
## 1 0.0000000 0.7777778 0.7500000 0.7142857
## 2 0.7777778 0.0000000 0.6000000 1.0000000
## 3 0.7500000 0.6000000 0.0000000 0.7500000
## 4 0.7142857 1.0000000 0.7500000 0.0000000

Remarquez que vegdist retourne une matrice dont la diagonale est de 0 (on l’affiche en spécifiant diag = TRUE). La diagonale est l’association d’un objet avec lui-même. Or la similarité d’un objet avec lui-même devrait être de 1! En fait, par convention vegdist retourne des dissimilarités, non pas des similarités. La matrice de distance serait donc calculée en extrayant la racine carrée des éléments de la matrice de dissimilarité:

##           1         2         3         4
## 1 0.0000000 0.8819171 0.8660254 0.8451543
## 2 0.8819171 0.0000000 0.7745967 1.0000000
## 3 0.8660254 0.7745967 0.0000000 0.8660254
## 4 0.8451543 1.0000000 0.8660254 0.0000000

Dans le chapitre sur l’analyse compositionnelle, nous avons abordé les significations différentes que peuvent prendre le zéro. L’information fournie par un zéro peut être différente selon les circonstances. Dans le cas d’une variable continue, un zéro signifie généralement une mesure sous le seuil de détection. Deux tissus dont la concentration en cuivre est nulle ont une afinité sous la perspective de la concentration en cuivre. Dans le cas de mesures d’abondance (décompte) ou d’occurence (présence-absence), on pourra décrire comme similaires deux niches écologiques où l’on retrouve une espèce en particulier. Mais deux sites où l’on de retouve pas d’ours polaires ne correspondent pas nécessairement à des niches similaires! En effet, il peut exister de nombreuses raisons écologiques et méthodologiques pour lesquelles l’espèces ou les espèces n’ont pas été observées. C’est le problème des double-zéros (espèces non observées à deux sites), problème qui est amplifié avec les grilles comprenant des espèces rares.

La ressemblance entre des objets comprenant des données continues devrait être calculée grâce à des indicateurs symétriques. Inversement, les affinités entre les objets décrits par des données d’abondance ou d’occurence susceptibles de générer des problèmes de double-zéros devraient être évaluées grâce à des indicateurs asymétriques. Un défi supplémentaire arrive lorsque les données sont de type mixte.

Nous utiliserons la convention de vegan et nous calculerons la dissimilarité, non pas la similarité. Les mesures de dissimilarité sont calculées sur des données d’abondance ou des données d’occurence. Notons qu’il existe beaucoup de confusion dans la littérature sur la manière de nommer les dissimilarités (ce qui n’est pas le cas des distances, dont les noms sont reconnus). Dans les sections suivantes, nous noterons la dissimilarité avec un \(d\) minuscule et la distance avec un \(D\) majuscule.

9.2.1 Association entre objets (mode Q)

9.2.1.1 Objets: Abondance

La dissimilarité de Bray-Curtis est asymétrique. Elle est aussi appelée l’indice de Steinhaus, de Czekanowski ou de Sørensen. Il est important de s’assurer de bien s’entendre la méthode à laquelle on fait référence. L’équation enlève toute ambiguité. La dissimilarité de Bray-Curtis entre les points A et B est calculée comme suit.

\[d_{AB} = \frac {\sum \left| A_{i} - B_{i} \right| }{\sum \left(A_{i}+B_{i}\right)}\]

Utilisons vegdist pour générer les matrices d’association. Le format “liste” de R est pratique pour enregistrer la collection d’objets, dont les matrice d’association que nous allons créer dans cette section.

##           1         2         3
## 2 0.9433962                    
## 3 0.9619048 0.4400000          
## 4 0.9591837 1.0000000 0.7647059

La dissimilarité de Bray-Curtis est souvent utilisée dans la littérature. Toutefois, la version originale de Bray-Curtis n’est pas tout à fait métrique (semimétrique). Conséquemment, la dissimilarité de Ruzicka (une variante de la dissimilarité de Jaccard pour les données d’abondance) est métrique, et devrait probablement être préféré à Bary-Curtis (Oksanen, 2006).

\[d_{AB, Ruzicka} = \frac { 2 \times d_{AB, Bray-Curtis} }{1 + d_{AB, Bray-Curtis}}\]

La dissimilarité de Kulczynski (aussi écrit Kulsinski) est asymétrique et semimétrique, tout comme celle de Bray-Curtis. Elle est calculée comme suit.

\[d_{AB} = 1-\frac{1}{2} \times \left[ \frac{\sum min(A_i, B_i)}{\sum A_i} + \frac{\sum min(A_i, B_i)}{\sum B_i} \right]\]

Une approche commune pour mesurer l’association entre sites décrits par des données d’abondance est la distance de Hellinger. Notez qu’il s’agit ici d’une distance, non pas d’une dissimilarité. Pour l’obtenir, on doit d’abord diviser chaque donnée d’abondance par l’abondance totale pour chaque site pour obtenir les espèces en tant que proportions, puis on extrait la racine carrée de chaque élément. Enfin, on calcule la distance euclidienne entre les proportions de chaque site. Pour rappel, une distance euclidienne est la généralisation en plusieurs dimensions du théorème de Pythagore, \(c = \sqrt{a^2 + b^2}\).

\[D_{AB} = \sqrt {\sum \left( \frac{A_i}{\sum A_i} - \frac{B_i}{\sum B_i} \right)^2}\]

😱 Attention La distance d’Hellinger hérite des biais liées aux données compositionnelles. Elle peut être substitiée par une matrice de distances d’Aitchison.

Toute comme la distance d’Hellinger, la distance de chord est calculée par une distance euclidienne sur des données d’abondance transformées de sorte que chaque ligne ait une longueur (norme) de 1.

La métrique du chi-carré, ou \(\chi\)-carré, ou chi-square, donne davantage de poids aux espèces rares qu’aux espèces communes. Son utilisation est recommandée lorsque les espèces rares sont de bons indicateurs de conditions écologiques particulières (Legendre et Legendre, 2012, p. 308).

\[ d_{AB} = \sqrt{\sum _j \frac{1}{\sum y_j} \left( \frac{A_j}{\sum A} - \frac{B_j}{\sum B} \right)^2 } \]

La métrique peut être transformée en distance en la multipliant par la racine carrée de la somme totale des espèces dans la matric d’abondance (\(X\)).

\[ D_{AB} = \sqrt{\sum X} \times d_{AB} \]

Une mannière visuellement plus intéressante de présenter une matrice d’association est un graphique de type heatmap.

Peu importe le type d’association utilisée, les heatmaps montrent les mêmes tendances. Les assocaitions de dissimilarité (Bray-Curtis, Kulczynski et Ruzicka) s’étalent de 0 à 1, tandis que les distances (Chi-Square, Chord et Hellinger) partent de zéro, mais n’ont pas de limite supérieure. On note les plus grandes différences entre les sites 2 et 4, tandis que les sites 2 et 3 sont les plus semblables pour toutes les mesures d’association à l’exception de la dissimilarité de Kulczynski.

9.2.1.2 Objets: Occurence (présence-absence)

Des indices d’association différents devraient être utilisés lorsque des données sont compilées sous forme booléenne. En général, les tableaux de données d’occurence seront compilés avec des 1 (présence) et des 0 (absence).

La similarité de Jaccard entre le site A et le site B est la proportion de double 1 (présences de 1 dans A et B) parmi les espèces. La dissimilarié est la proportion complémentaire (comprenant [1, 0], [0, 1] et [0, 0]). La distance de Jaccard est la racine carrée de la dissimilarité.

Les distances d’Hellinger, de chord et de chi-carré sont aussi appropriées pour les calculs de distances sur des tableaux d’occurence.

Graphiquement,

Il est attendu que les matrices d’association sur l’occurence sont semblables à celles sur l’abondance. Dans ce cas-ci, la distance d’Hellinger donne des résultats semblables à la dissimilarité de Jaccard.

9.2.1.3 Objets: Données quantitatives

Les données quantitative en écologie peuvent décrire l’état de l’environnement: le climat, l’hydrologie, l’hydrogéochimie, la pédologie, etc. En règle générale, les coordonnées des sites ne sot pas des variables environnementales, à que l’on soupçonne la coordonnée elle-même d’être responsable d’effets sur notre système: mais il s’agira la plupart du temps d’effets confondants (par exemple, on peut mesurer un effet de lattitude sur le rendement des agrumes, mais il s’agira probablement avant tout d’effets dus aux conditions climatiques, qui elles changent en fonction de la lattitude). D’autre types de données quantitative pouvant être appréhendées par des distances sont les traits phénologiques, les ionomes, les génomes, etc.

La distance euclidienne est la racine carrée de la somme des carrés des distances sur tous les axes. Il s’agit d’une application multidimensionnelle du théorème de Pythagore. La distance d’Aitchison, couverte dans le chapitre 8, est une distance euclidienne calculée sur des données compositionnelles préalablement transformées. La distance euclidienne est sensible aux unités utilisés: utiliser des milimètres plutôt que des mètres enflera la distance euclidienne. Il est recommandé de porter une attention particulière aux unités, et de standardiser les données au besoin (par exemple, en centrant la moyenne à zéro et en fixant l’écart-type à 1).

On pourrait, par exemple, mesurer la distance entre des observations des dimensions de différentes espèces d’iris. Ce tableau est inclu dans R par défaut.

##   Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
## 1          6.6         3.0          4.4         1.4 versicolor
## 2          7.7         3.8          6.7         2.2  virginica
## 3          4.4         3.2          1.3         0.2     setosa
## 4          5.6         3.0          4.1         1.3 versicolor
## 5          6.9         3.1          5.1         2.3  virginica

Les mesures du tableau sont en centimètres. Pour éviter de donner davantage de poids aux longueur des sépales et en même temps de négliger la largeur des pétales, nous allons standardiser le tableau.

## # A tibble: 150 x 5
##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
##           <dbl>       <dbl>        <dbl>       <dbl> <fct>  
##  1       -0.898      1.02          -1.34       -1.31 setosa 
##  2       -1.14      -0.132         -1.34       -1.31 setosa 
##  3       -1.38       0.327         -1.39       -1.31 setosa 
##  4       -1.50       0.0979        -1.28       -1.31 setosa 
##  5       -1.02       1.25          -1.34       -1.31 setosa 
##  6       -0.535      1.93          -1.17       -1.05 setosa 
##  7       -1.50       0.786         -1.34       -1.18 setosa 
##  8       -1.02       0.786         -1.28       -1.31 setosa 
##  9       -1.74      -0.361         -1.34       -1.31 setosa 
## 10       -1.14       0.0979        -1.28       -1.44 setosa 
## # ... with 140 more rows

Pour les comparaisons des dimensions, prenons la moyenne des dimensions (mises à l’échelle) par espèce.

## # A tibble: 3 x 4
##   Sepal.Length Sepal.Width Petal.Length Petal.Width
##          <dbl>       <dbl>        <dbl>       <dbl>
## 1       -1.01        0.850       -1.30       -1.25 
## 2        0.112      -0.659        0.284       0.166
## 3        0.899      -0.191        1.02        1.08

Nous pouvons utiliser la distance euclidienne, commune en géométrie, pour comparer les espèces. La distance euclidienne est calculée comme suit.

\[ \mathcal{E} = \sqrt{\Sigma_i \left( A_i - B_i \right) ^2 } \]

La distance de Mahalanobis est semblable à la distance euclidienne, mais qui tient compte de la covariance de la matrice des objets. Cette covariance peut être utilisée pour décrire la structure d’un nuage de points. La diastance de Mahalanobis se calcule comme suit.

\[\mathcal{M} = \sqrt{(A - B)^T S^{-1} (A-B)}\]

Notez qu’il s’agit d’une généralisation de la distance euclidienne, qui équivaut à une distance de Mahalanobis dont la matrice de covariance est une matrice identité.

La distance de Mahalanobis permet de représenter des distances dans un espace fortement corrélé. Elle est courramment utilisée pour détecter les valeurs aberrantes selon des critères de distance à partir du centre d’un jeu de données multivariées.

La distance de Manhattan porte aussi le nom de distance de cityblock ou de taxi. C’est la distance que vous devrez parcourir pour vous rendre du point A au point B à Manhattan, c’est-à-dire selon une séquence de tronçons perpendiculaires.

\[ D_{AB} = \sum _i \left| A_i - B_i \right| \]

La distance de Manhattan est appropriée lorsque les gradients (changements d’un état à l’autre ou d’une région à l’autre) ne permettent pas des changements simultanés. Mieux vaut standardiser les variables pour éviter qu’une dimension soit prépondérante.

Avant de présenter les résultats des espèces d’iris, voici une représentation des distances euclidiennes (rouge), de Mahalanobis (bleu) et de Manhattan (vert), chacune de 1 et 2 unités à partir du centre et, pour ce qui est de la distance de Mahalanobis, selon la covariance.

## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select

Et, graphiquement, les résultats des distances des iris.

Le tableau iris est ordonné par espèce. Les distances euclidienne et de Manhattan permettent aisément de distinguer les espèces selon les dimensions des pétales et des sépales. Toutefois, l’utilsation de la covariance avec la distance de Mahalanobis crée des distinction moins tranchées.

9.2.1.4 Objets: Données mixtes

Les données catégorielles ordinales peuvent être transformées en données continues par gradations linéaires ou quadratiques. Les données catégorielles nominales, quant à elles, peuvent être encodées (encodage catégoriel) en données similaires à des occurences. Attention toutefois: contrairement à la régression linéaire qui demande d’exclure une catégorie, l’encodage catégoriel doit inclure toutes les catégories. Le comportement par défaut de la fonction model.matrix est d’exclure la catégorie de référence: on doit spécifier que l’intercept est de zéro, c’est-à-dire model.matrix(~ + categorie).

La similarité de Gower a été développée pour mesurer des associations entre des objets dont les données sont mixtes: booléennes, catégorielles et continues. La similarité de Gower est calculée en additionnant les distances calculées par colonne, individuellement. Si la colonne est booléenne, on utilise les distances de Jaccard (qui exclue les double-zéro) de manière univariée: une variable à la fois. Pour les variables continues, on utilise la distance de Manhattan divisée par la plage de valeurs de la variable (pour fin de standardisation). Puisqu’elle hérite de la particularité de la distance de Manhattan et de la similarité de Jaccard univariée, la similarité de Gower reste une combinaison linéaire de distances univariées.

## # A tibble: 8 x 7
##      ID   age gender civil_status salary children available_credit
##   <int> <dbl> <chr>  <chr>         <dbl> <lgl>               <dbl>
## 1     1    21 M      MARRIED        3000 TRUE                 2200
## 2     2    21 M      SINGLE         1200 FALSE                 100
## 3     3    19 N      SINGLE        32000 TRUE                22000
## 4     4    30 M      SINGLE         1800 TRUE                 1100
## 5     5    21 F      MARRIED        2900 TRUE                 2000
## 6     6    21 F      SINGLE         1100 TRUE                  100
## 7     7    19 F      WIDOW         10000 FALSE                6000
## 8     8    30 F      DIVORCED       1500 TRUE                 2200

Il faut préalablement procéder à l’encodage catégoriel pour les variables catégorielles nominales.

##   age genderF genderM genderN civil_statusMARRIED civil_statusSINGLE civil_statusWIDOW salary
## 1  21       0       1       0                   1                  0                 0   3000
## 2  21       0       1       0                   0                  1                 0   1200
## 3  19       0       0       1                   0                  1                 0  32000
## 4  30       0       1       0                   0                  1                 0   1800
## 5  21       1       0       0                   1                  0                 0   2900
## 6  21       1       0       0                   0                  1                 0   1100
## 7  19       1       0       0                   0                  0                 1  10000
## 8  30       1       0       0                   0                  0                 0   1500
##   childrenTRUE available_credit
## 1            1             2200
## 2            0              100
## 3            1            22000
## 4            1             1100
## 5            1             2000
## 6            1              100
## 7            0             6000
## 8            1             2200
## attr(,"assign")
##  [1] 1 2 2 2 3 3 3 4 5 6
## attr(,"contrasts")
## attr(,"contrasts")$gender
## [1] "contr.treatment"
## 
## attr(,"contrasts")$civil_status
## [1] "contr.treatment"
## 
## attr(,"contrasts")$children
## [1] "contr.treatment"

Calculons la dissimilarité de Gower (cette fois le graphique est fait avec pheatmap).

Les dendrogrammes apparaissants sur les axes du graphique sont issus d’un processus de partitionnement basé sur la distance, que nous verrons plus loin dans ce chapiter. Les profils des clients 4 et 7, ainsi que ceux des clients 3 et 7 diffèrent le plus. Les profils 3 et 4 sont néanmoins plutôt différents.

9.2.2 Associations entre variables (mode R)

Il existe de nombreuses approches pour mesurer les associations entre variables. La plus connue est la corrélation. Mais les données d’abondance et d’occurence demandent des approches différentes.

9.2.2.1 Variables: Abondance

La distance du chi-carré est suggérée par Borcard et al. (2011).

Des coabondances sont notables pour la mésange à tête noire, le jaseur boréal, la citelle à poitrine rousse et le bruant à gorge blanche (tache bleu au centre).

9.2.2.2 Variables: Occurence

La dissimilarité de Jaccard peut être utilisée.

Des cooccurences sont notables pour le jaseur boréal, la citelle à poitrine rousse et le bruant à gorge blanche (tache bleu au centre).

9.2.2.3 Variables: Quantités

La matrice des corrélations de Pearson peut être utilisée pour les données continues. Quant aux variables ordinales, elles devraient idéalement être liées linéairement ou quadratiquement. Si ce n’est pas le cas, c’est-à-dire que les catégories sont ordonnées par rang seulement, vous pourrez avoir recours aux coefficients de corrélation de Spearman ou de Kendall.

9.2.3 Conclusion sur les associations

Il n’existe pas de règle claire pour déterminer quelle technique d’association utiliser. Cela dépend en premier lieu de vos données. Vous sélectionnerez votre méthode d’association selon le type de données que vous abordez, la question à laquelle vous désirez répondre ainsi l’expérience dans la littérature comme celle de vos collègues scientifiques. S’il n’existe pas de règle clair, c’est qu’il existe des dizaines de méthodes différentes, et la plupart d’entre elles vous donneront une perspective juste et valide. Il faut néanmoins faire attention pour éviter de sélectionner les méthodes qui ne sont pas appropriées.

9.3 Partitionnement

Les données suivantes ont été générées par Leland McInnes (Tutte institute of mathematics, Ottawa). Êtes-vous en mesure d’identifier des groupes? Combien en trouvez-vous?

## Parsed with column specification:
## cols(
##   x = col_double(),
##   y = col_double()
## )

En 2D, l’oeil humain peut facilement détecter les groupes. En 3D, c’est toujours possible, mais au-delà de 3D, le partitionnement cognitive devient rapidement maladroite. Les algorithmes sont alors d’une aide précieuse. Mais ils transportent en pratique tout un baggage de limitations. Quel est le critère d’association entre les groupes? Combien de groupe devrions-nous créer? Comment distinguer une donnée trop bruitée pour être classifiée?

Le partitionnement de données (clustering en anglais), et inversement leur regroupement, permet de créer des ensembles selon des critères d’association. On suppose donc que Le partitionnement permet de créer des groupes selon l’information que l’on fait émerger des données. Il est conséquemment entendu que les données ne sont pas catégorisées à priori: il ne s’agit pas de prédire la catégorie d’un objet, mais bien de créer des catégories à partir des objets par exemple selon leurs dimensions, leurs couleurs, leurs signature chimique, leurs comportements, leurs gènes, etc.

Plusieurs méthodes sont aujourd’hui offertes aux analystes pour partitionner leurs données. Dans le cadre de ce manuel, nous couvrirons ici deux grandes tendances dans les algorithmes.

  1. Méthodes hiérarchique et non hiérarchiques. Dans un partitionnement hiérarchique, l’ensemble des objets forme un groupe, comprenant des sous-regroupements, des sous-sous-regroupements, etc., dont les objets forment l’ultime partitionnement. On pourra alors identifier comment se décline un partitionnement. À l’inverse, un partitionnement non-hiérarchique des algorhitmes permettent de créer les groupes non hiérarchisés les plus différents que possible.

  2. Membership exclusif ou flou. Certaines techniques attribuent à chaque objet une classe unique: l’appartenance sera indiquée par un 1 et la non appartenance par un 0. D’autres techniques vont attribuer un membership flou où le degré d’appartenance est une variable continue de 0 à 1. Parmi les méthodes floues, on retrouve les méthodes probabilistes.

9.3.1 Évaluation d’un partitionnement

Le choix d’une technique de partitionnement parmi de nombreuses disponibles, ainsi que le choix des paramètres gouvernant chacune d’entre elles, est avant tout basé sur ce que l’on désire définir comme étant un groupe, ainsi que la manière d’interpréter les groupes. En outre, le nombre de groupe à départager est toujours une décision de l’analyste. Néanmoins, on peut se fier des indicateurs de performance de partitionnement. Parmis ceux-ci, retenons le score silouhette ainsi que l’indice de Calinski-Harabaz.

9.3.1.1 Score silouhette

En anglais, le h dans silouhette se trouve après le l: on parle donc de silhouette coefficient pour désigner le score de chacun des objets dans le partitionnement. Pour chaque objet, on calcule la distance moyenne qui le sépare des autres points de son groupe (\(a\)) ainsi que la distance moyenne qui le sépare des points du groupe le plus rapproché.

\[s = \frac{b-a}{max \left(a, b \right)}\]

Un coefficient de -1 indique le pire classement, tandis qu’un coefficient de 1 indique le meilleur classement. La moyenne des coefficients silouhette est le score silouhette.

9.3.1.2 Indice de Calinski-Harabaz

L’indice de Calinski-Harabaz est proportionnel au ratio des dispersions intra-groupe et la moyenne des dispersions inter-groupes. Plus l’indice est élevé, mieux les groupes sont définis. La mathématique est décrite dans la documentation de scikit-learn, un module d’analyse et autoapprentissage sur Python.

Note. Les coefficients silouhette et l’indice de Calinski-Harabaz sont plus appropriés pour les formes de groupes convexes (cercles, sphères, hypersphères) que pour les formes irrégulières (notamment celles obtenues par la DBSCAN, discutée ci-desssous).

9.3.2 Partitionnement non hiérarchique

Il peut arriver que vous n’ayez pas besoin de comprendre la structure d’agglomération des objets (ou variables). Plusieurs techniques de partitionnement non hiérarchique sont disponibles sur R. On s’intéressera en particulier aux k-means et au dbscan.

9.3.2.1 Kmeans

L’objectif des kmeans est de minimiser la distance euclédienne entre un nombre prédéfini de k groupes exclusifs.

  1. L’algorhitme commence par placer une nombre k de centroides au hasard dans l’espace d’un nombre p de variables (vous devez fixer k, et p est le nombre de colonnes de vos données).
  2. Ensuite, chaque objet est étiquetté comme appartenant au groupe du centroid le plus près.
  3. La position du centroide est déplacée à la moyenne de chaque groupe.
  4. Recommencer à partir de l’étape 2 jusqu’à ce que l’assignation des objets aux groupes ne change plus.
Source: David Sheehan

La technique des kmeans suppose que les groupes ont des distributions multinormales - représentées par des cercles en 2D, des sphères en 3D, des hypersphères en plus de 3D. Cette limitation est problématique lorsque les groupes se présentent sous des formes irrégulières, comme celles du nuage de points de Leland McInnes, présenté plus haut. De plus, la technique classique des kmeans est basée sur des distances euclidiennes: l’utilisation des kmeans n’est appropriée pour les données comprenant beaucoup de zéros, comme les données d’abondance, qui devraient préalablement être transformées en variables centrées et réduites (Legendre et Legendre, 2012). La technique des mixtures gaussiennes (gaussian mixtures) est une généralisation des kmeans permettant d’intégrer la covariance des groupes. Les groupes ne sont plus des hyper-sphères, mais des hyper-ellipsoïdes.

9.3.2.1.1 Application

Nous pouvons utilisé la fonction kmeans de R. Toutefois, puisque l’on désire ici effectuer des tests de partitionnement pour plusieurs nombres de groupes, nous utiliserons cascadeKM, du module vegan. Notez que de nombreux paramètres par défaut sont utilisés dans les exécutions ci-dessous. Ces notes de cours ne forment pas un travail de recherche scientifique. Lors de travaux de recherche, l’utilsation d’un argument ou d’un autre dans une fonction doit être justifié: qu’un paramètre soit utilisé par défaut dans une fonction n’est a priori pas une justification convainquante.

Pour les kmeans, on doit fixer le nombre de groupes. Le graphique des données de Leland McInnes montrent 6 groupes. Toutefois, il est rare que l’on puisse visualiser des démarquations aussi tranchées que celles de l’exemple, qui plus est dans des cas où l’on doit traiter de plus de deux dimensions. Je vais donc lancer le partitionnement en boucle pour plusieurs nombres de groupes, de 3 à 10 et pour chaque groupe, évaluer le score silouhette et de Calinski-Habaraz. J’utilise un argument random_state pour m’assurer que les groupes seront les mêmes à chaque fois que la cellule sera lancée.

## List of 4
##  $ partition: int [1:2309, 1:8] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2309] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:8] "3 groups" "4 groups" "5 groups" "6 groups" ...
##  $ results  : num [1:2, 1:8] 85.1 2164.5 61.4 2294.6 51.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "SSE" "calinski"
##   .. ..$ : chr [1:8] "3 groups" "4 groups" "5 groups" "6 groups" ...
##  $ criterion: chr "calinski"
##  $ size     : int [1:10, 1:8] 1243 505 561 NA NA NA NA NA NA NA ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:10] "Group 1" "Group 2" "Group 3" "Group 4" ...
##   .. ..$ : chr [1:8] "3 groups" "4 groups" "5 groups" "6 groups" ...
##  - attr(*, "class")= chr "cascadeKM"

L’objet mcinnes_kmeans, de type cascadeKM, peut être visualisé directement avec la fonction plot.

On obtient un maximum de Calinski à 4 groupes, qui correspons à la deuxième simulation effectuée de 3 à 10.

Examinons les scores silouhette (module: cluster).

Le score silouhette maximum est à 3 groupes. La forme des groupes n’étant pas convexe, il fallait s’attendre à ce que indicateurs maximaux pour les deux indicateurs soient différents. C’est d’ailleurs souvent le cas. Cet exemple supporte que le choix du nombre de groupe à départager repose sur l’analyste, non pas uniquement sur les indicateurs de performance. Choisissons 6 groupes, puisque que c’est visuellement ce que l’on devrait chercher pour ce cas d’étude.

##   3 groups 4 groups 5 groups 6 groups 7 groups 8 groups 9 groups 10 groups
## 1        1        4        5        3        1        6        7        10
## 2        1        4        4        1        1        3        7        10
## 3        1        1        5        3        6        6        3         3

L’algorithme kmeans est loin d’être statisfaisant. Cela est attendu, puisque les kmeans recherchent des distribution gaussiennes sur des groupes vraisemblablement non-gaussiens.

Nous pouvons créer un graphique silouhette pour nos 6 groupes. Notez qu’à cause d’un bogue, il n’est pas possible de présenter les données clairement lorsqu’elles sont nombreuses.

9.3.2.2 DBSCAN

La technique DBSCAN (* Density-Based Spatial Clustering of Applications with Noise) sousentend que les groupes sont composés de zones où l’on retrouve plus de points (zones denses) séparées par des zones de faible densité. Pour lancer l’algorithme, nous devons spécifier une mesure d’association critique (distance ou dissimilarité) d* ainsi qu’un nombre de point critique k dans le voisinage de cette distance.

  1. L’algorithme commence par étiqueter chaque point selon l’une de ces catégories:
  • Noyau: le point a au moins k points dans son voisinage, c’est-à-dire à une distance inférieure ou égale à d.
  • Bordure: le point a moins de k points dans son voisinage, mais l’un de des points voisins est un noyau.
  • Bruit: le cas échéant. Ces points sont considérés comme des outliers.

  1. Les noyaux distancés de d ou moins sont connectés entre eux en englobant les bordures.

Le nombre de groupes est prescrit par l’algorithme DBSCAN, qui permet du coup de détecter des données trop bruitées pour être classées.

Damiani et al. (2014) a développé une approche utilisant la technique DBSCAN pour partitionner des zones d’escale pour les flux de populations migratoires.

9.3.2.2.1 Application

La technique DBSCAN n’est pas basée sur le nombre de groupe, mais sur la densité des points. L’argument x ne constitue pas les données, mais une matrice d’association. L’argument minPts spécifie le nombre minimal de points qui l’on doit retrouver à une distance critique d* pour la formation des *noyaux et la propagation des groupes, spécifiée dans l’argument eps. La distance d peut être estimée en prenant une fraction de la moyenne, mais on aura volontiers recours à sont bon jugement.

## [1] 1 0 2 6 3 4 5

Les paramètres spécifiés donnent 5 groupes (1, 2, ..., 5) et des points trop bruités pour être classifiés (étiquetés 0). Voyons comment les groupes ont été formés.

Le partitionnement semble plus conforme à ce que l’on recherche. Néanmoins, DBSCAN cré quelques petits groupes indésirables (groupe 6, en rose) ainsi qu’un grand groupe (violet) qui auraient lieu d’être partitionné. Ces défaut pourraient être réglés en jouant sur les paramètres eps et minPts.

9.3.3 Partitionnement hiérarchique

Les techniques de partitionnement hiérarchique sont basées sur les matrices d’association. La technique pour mesurer l’association (entre objets ou variables) déterminera en grande partie le paritionnement des données. Les partitionnements hiérarchiques ont l’avantage de pouvoir être représentés sous forme de dendrogramme (ou arbre) de partition. Un tel dendrogramme présente des sous-groupes qui se joignent en groupes jusqu’à former un seul ensemble.

Le partitionnement hiérarchique est abondamment utilisé en phylogénie, pour étudier les relations de parenté entre organismes vivants, populations d’organismes et espèces. La phénétique, branche empirique de la phylogénèse interspécifique, fait usage du partitionnement hiérarchique à partir d’associations génétiques entre unités taxonomiques. On retrouve de nombreuses ressources académiques en phylogénétique ainsi que des outils pour R et Python. Toutefois, la phylogénétique en particulier ne fait pas partie de la présente ittération de ce manuel.

9.3.3.1 Techniques de partitionnement hiérarchique

Le partitionnement hiérarchique est typiquement effectué avec une des quatres méthodes suivantes, dont chacune possède ses particularités, mais sont toutes agglomératives: à chaque étape d’agglomération, on fusionne les deux groupes ayant le plus d’affinité sur la base des deux sous-groupes les plus rapprochés.

Single link (single). Les groupes sont agglomérés sur la base des deux points parmi les groupes, qui sont les plus proches.

Complete link (complete). À la différence de la méthode single, on considère comme critère d’agglomération les éléments les plus éloignés de chaque groupe.

Agglomération centrale. Il s’agit d’une fammlle de méthode basées sur les différences entre les tendances centrales des objets ou des groupes.

  • Average (average). Appelée UPGMA (Unweighted Pair-Group Method unsing Average), les groupes sont agglomérés selon un centre calculés par la moyenne et le nombre d’objet pondère l’agglomération (le poids des groupes est retiré). Cette technique est historiquement utilisée en bioinformatique pour partitionner des groupes phylogénétiques (Sneath et Sokal, 1973).
  • Weighted (weighted). La version de average, mais non pondérée (WPGMA).
  • Centroid (centroid). Tout comme average, mais le centroïde (centre géométrique) est utilisé au lieu de la moyenne. Accronyme: UPGMC.
  • Median (median). Appelée WPGMC. Devinez! ;)

Ward (ward). L’optimisation vise à minimiser les sommes des carrés par regroupement.

9.3.3.2 Quel outil de partitionnement hiérarchique utiliser?

Alors que le choix de la matrice d’association dépend des données et de leur contexte, la technique de partitionnement hiérarchique peut, quant à elle, être basée sur un critère numérique. Il en existe plusieurs, mais le critère recommandé pour le choix d’une technique de partitionnement hiérarchique est la corrélation cophénétique. La distance cophénétique est la distance à laquelle deux objets ou deux sous-groupes deviennent membres d’un même groupe. La corrélation cophénétique est la corrélation de Pearson entre le vecteur d’association des objets et le vecteur de distances cophénétiques.

9.3.3.3 Application

Les techniques de partitionnement hiérarchique présentées ci-dessus sont disponibles dans le module stats de R, qui est chargé automatiquement lors de l’ouversture de R. Nous allons classifier les dimensions des iris grâce à la distance de Manhattan.

## The "ward" method has been renamed to "ward.D"; note new "ward.D2"

La méthode average retourne la corrélation la plus élevée. Pour plus de flexibilité, enchâssons le nom de la méthode dans une variable. Ainsi, en chageant le nom de cette variable, le reste du code sera conséquent.

Le partitionnement hiérarchique peut être visualisé par un dendrogramme.

9.3.3.4 Combien de groupes utiliser?

La longueur des lignes verticales est la distance séparant les groupes enfants. Bien que la sélection du nombre de groupe soit avant tout basée sur les besoins du problème, nous pouvons nous appuyer sur certains outils. La hauteur totale peut servir de critère pour définir un nombre de groupes adéquat. On pourra sélectionner le nombre de groupe où la hauteur se stabilise en fonction du nombre de groupe. On pourra aussi utiliser le graphique silhouette, comprenant une collection de largeurs de silouhette, représentant le degré d’appartenance à son groupe. La fonction sklearn.metrics.silhouette_score, du module scikit-learn, s’en occupe.

Le nombre optimal de groupes serait de 5. Coupons le dendrorgamme à la hauteur correspondant à 5 groupes avec la fonction cutree.

La classification hiérarchique, uniquement basée sur la distance, peut être inappropriée pour définir des formes complexes.

9.3.4 Partitionnement hiérarchique basée sur la densité des points

La tecchinque HDBSCAN, dont l’algorithme est relativement récent (Campello et al., 2013), permet une partitionnement hiérarchique sur le même principe des zones de densité de la technique DBSCAN. Le HDBSCAN a été utilisée pour partitionner les lieux d’escale d’oiseaux migrateurs en Chine (Xu et al., 2013).

Avec DBSCAN, un rayon est fixé dans une métrique appropriée. Pour chaque point, on compte le nombre de point voisins, c’est à dire le nombre de point se situant à une distance (ou une dissimilarité) égale ou inférieure au rayon fixé. Avec HDBSCAN, on spécifie le nombre de points devant être recouverts et on calcule le rayon nécessaire pour les recouvrir. Ainsi, chaque point est associé à un rayon critique que l’on nommera \(d_{noyau}\). La métrique initiale est ensuite altérée: on remplace les associations entre deux objets A et B par la valeur maximale entre cette association, le rayon critique de A et le rayon critique de B. Cette nouvelle distance est appelée la distance d’atteinte mutuelle: elle accentue les distances pour les points se trouvant dans des zones peu denses. On applique par la suite un algorithme semblable à la partition hiérarchique single link: En s’élargissant, les rayons se superposent, chaque superposition de rayon forment graduellement des groupes qui s’agglomèrent ainsi de manière hiérarchique. Au lieu d’effectuer une tranche à une hauteur donnée dans un dendrogramme de partitionnement, la technique HDBSCAN se base sur un dendrogramme condensé qui discarte les sous-groupes comprenant moins de n objets (\(n_{gr min}\)). Dans nouveau dendrogramme, on recherche des groupes qui occupent bien l’espace d’analyse. Pour ce faitre, on utilise l’inverse de la distance pour créer un indicateur de persistance (semblable à la similarité), \(\lambda\). Pour chaque groupe hiérarchique dans le dendrogramme condensé, on peut calculer la persistance où le groupe prend naissance. De plus, pour chaque objet d’un groupe, on peut aussi calculer une distance à laquelle il quitte le groupe. La stabilité d’un groupe est la domme des différences de persistance entre la persistance à la naissance et les persistances des objets. On descend dans le dendrogramme. Si la somme des stabilité des groupes enfants est plus grande que la stabilité du groupe parent, on accepte la division. Sinon, le parent forme le groupe. La documentation du module hdbscan pour Python offre une description intuitive et plus exhaustive des principes et algorithme de HDBSCAN.

9.3.4.1 Paramètres

Outre la métrique d’association dont nous avons discuté, HDBSCAN demande d’être nourri avec quelques paramètres importants. En particulier, le nombre minimum d’objets par groupe, \(n_{gr min}\) dépend de la quantité de données que vous avez à votre disposition, ainsi que de la quantité d’objets que vous jugez suffisante pour créer des groupes. Nous utiliserons l’implémentation de HDBSCAN du module dbscan. Si vous désirez davantage d’options, vous préférerez probablement l’implémentation du module largeVis.

## [1] 6 0 4 3 5 1 2

Nous avons 6 groupes, numérotés de 1 à 6, ainsi que des étiquettes identifiant des objets désignés comme étant du bruit de fond, numéroté 0. Le dendrogramme non condensé peu être produit.

Difficile d’y voir clair avec autant d’objets. L’objet mcinnes_hdbscan a un nombre minimum d’objets par groupe de 20. Ce qui permet de présenter le dendrogramme de manière condensée.

Enfin, un aperçu des stratégies de partitionnement utilisés jusqu’ici.

## Warning: attributes are not identical across measure variables;
## they will be dropped

Clairement, le partitionnement avec HDBSCAN donne les meilleurs résultats.

9.3.5 Conclusion sur le partitionnement

Au chapitre 4, nous avons vu avec le jeu de données “datasaurus” que la visualisation peut permettre de détecter des structures en segmentant les données selon des groupes.

Or, si les données n’étaient pas étiquetées, leur structure serait indétectable avec les algorithmes disponibles actuellement. Le partitionnement permet d’explorer des données, de détecter des tendances et de dégager des groupes permettant la prise de décision.

Plusieurs techniques de partitionnement ont été présentées. Le choix de la technique sera déterminante sur la manière dont les groupes seront partitionnés. La définition d’un groupe variant d’un cas à l’autre, il n’existe pas de règle pour prescrire une méthode ou une autre. La partitionnement hiérarchique a l’avantage de permetre de visualiser comment les groupes s’agglomèrent. Parmi les méthodes de partitionnement hiérarchique disponibles, les méthodes basées sur la densité permettent une grande flexibilité, ainsi qu’une détection d’observations ne faisant partie d’aucun goupe.

9.4 Ordination

En écologie, biologie, agronommie comme en foresterie, la plupart des tableaux de données comprennent de nombreuses variables: pH, nutriments, climat, espèces ou cultivars, etc. L’ordination vise à mettre de l’ordre dans des données dont le nombre élevé de variables peut amener à des difficultés d’appréciation et d’interprétaion (Legendre et Legendre, 2012). Plus précisément, le terme ordination est utilisé en écologie pour désigner les techniques de réduction d’axe. L’analyse en composante principale est probablement la plus connue de ces techniques. Mais de nombreuses techniques d’ordination ont été développées au cours des dernières années, chacune ayant ses domaines d’application.

Les techniques de réduction d’axe permettent de dégager l’information la plus importante en projetant une synthèse des relations entre les observations et entre les variables. Les techniques ne supposant aucune structure a priori sont dites non-contraignantes: elles ne comprennent pas de tests statistiques. À l’inverse, les ordinations contraignantes lient des variables descriptives avec une ou plusieurs variables prédictives.

La référence en la matière est indiscutablement (Legendre et Legendre, 2012). Cette section en couvrira quelques unes et vous guidera vers la technique la plus appropriée pour vos données.

9.4.1 Ordination non contraignante

Cette section couvrira l’analyse en composantes principales (ACP), l’analyse de correspondance (AC), l’analyse factorielle (AF) ainsi que l’analyse en coordonnées principales (ACoP).

Méthode Distance préservée Variables
Analyse en composantes principales (ACP) Distance euclidienne Données quantitatives, relations linéaires (attention aux double-zéros)
Analyse de correspondance (AC) Distance de \(\chi^2\) Données non-négatives, dimentionnellement homogènes ou binaires, abondance ou occurence
Positionnement multidimensionnel (PoMd) Toute mesure de dissimilarité Données quantitatives, qualitatives nominales/ordinales ou mixtes

Source: Adapté de (Legendre et Legendre, 2012, chapitre 9)

9.4.1.1 Analyse en composantes principales

L’objectif d’une ACP est de représenter les données dans un nombre réduit de dimensions représentant le plus possible la variation d’un tableau de données: elle permet de projetter les données dans un espace où les variables sont combinées en axes orthogonaux dont le premier axe capte le maximum de variance. L’ACP peut par exemple être utilisée pour analyser des corrélations entre variables ou dégager l’information la plus pertinente d’un tableau de données météo ou de signal en un nombre plus retreint de variables.

L’ACP effectue une rotation des axes à partir du centre (moyenne) du nuage de points effectuée de manière à ce que le premier axe définisse la direction où l’on retrouve la variance maximale. Ce premier axe est une combinaison linéaire des variables et forme la première composante principale. Une fois cet axe définit, on trouve de deuxième axe, orthogonal au premier, où l’on retouve la variance maximale - cet axe forme la deuxième composante principale, et ainsi de suite jusqu’à ce que le nombre d’axe corresponde au nombre de variables. Les projections des observations sur ces axes principaux sont appelés les scores. Les projections des variables sur les axes principaux sont les vecteurs propres (eigenvectors, ou loadings). La variance des composantes principales diminue de la première à la dernière, et peut être calculée comme une proportion de la variance totale: c’est le pourcentage d’inertie. Par convention, on utilise les valeurs propres (eigenvalues) pour mesurer l’importance des axes. Si la première composante principale a une inertie de 50% et la deuxième a une intertie de 30%, la représentation en 2D des projection représentera 80% de la variance du nuage de points.

L’hétérogénéité des échelles de mesure peut avoir une grande importance sur les résultats d’une ACP (les données doivent être dimensionnellement homogènes). En effet, la hauteur d’un ceriser aura une variance plus grande que le diamètre d’une cerise exprimé dans les mêmes unités, et cette dernière aura plus de variance que la teneur en cuivre d’une feuille. Il est conséquemment avisé de mettre les données à l’échelle en centrant la moyenne à zéro et l’écart-type à 1 avant de procéder à une ACP.

L’ACP a été conçue pour projetter en un nombre moindre de dimensions des observations dont les distributions sont multinormales. Bien que l’ACP soit une technique robuste, il est préférable de transformer préalablement les variables dont la distribution est particulièrement asymétriques (Legendre et Legendre, 2012, p. 450). Le cas échéant, les valeurs extrêmes pourraient faire dévier les vecteurs propres et biaiser l’analyse. En particulier, les données ACP menées sur des données compositionnelles sont réputées pour générer des analyses biaisées (Pawlowsky-Glahn and Egozcue, 2006). Le test de Mardia (Korkmaz, 2014) peut être utilisé pour tester la multinormalité. Une distribution multinormale devrait générer des scores en forme d’hypersphère (en forme de cercle sur un biplot: voir plus loin).

9.4.1.1.1 Vecteurs propres et valeurs propres

Une matrice carrée (comme une matrice de covariance \(\Sigma\)) multipliée par un vecteur propre \(e\) est égale aux valeurs propres \(\lambda\) multipliées par les vecteurs propres \(e\).

\[ \Sigma e = \lambda e \]

De manière intuitive, les vecteurs propres indiquent l’orientation de la covariance, et les valeurs propres indique la longueur associée à cette direction. L’ACP est basée sur le calcul des vecteurs propres et des valeurs propres de la matrice de covariance des variables. Pour d’abord obtenir les valeurs propres \(\lambda\), il faut résoudre l’équation

\[ det(cov(X) - \lambda I) = 0 \],

où \(det\) est l’opération permettant de calculer le déterminant, \(cov\) est l’opération pour calculer la covariance, \(X\) est la matrice de données, \(\lambda\) sont les valeurs propres et \(I\) est une matrice d’identité.

Pour \(p\) variables dans votre tableau \(X\), vous obtiendrex \(p\) valeurs propres. Ensuite, on trouve les vecteurs propres en résolvant l’équation $ e = e $.

Bien qu’il soit possible d’effectuer cette opération à la main pour des cas très simples, vous aurez avantage à utiliser un langage de programmation.

Chargeons les données d’iris, puis isolons seulement les deux dimensions des sépales l’espèce setosa.

##    Sepal.Length Sepal.Width
## 1           5.1         3.5
## 2           4.9         3.0
## 3           4.7         3.2
## 4           4.6         3.1
## 5           5.0         3.6
## 6           5.4         3.9
## 7           4.6         3.4
## 8           5.0         3.4
## 9           4.4         2.9
## 10          4.9         3.1
## 11          5.4         3.7
## 12          4.8         3.4
## 13          4.8         3.0
## 14          4.3         3.0
## 15          5.8         4.0
## 16          5.7         4.4
## 17          5.4         3.9
## 18          5.1         3.5
## 19          5.7         3.8
## 20          5.1         3.8
## 21          5.4         3.4
## 22          5.1         3.7
## 23          4.6         3.6
## 24          5.1         3.3
## 25          4.8         3.4
## 26          5.0         3.0
## 27          5.0         3.4
## 28          5.2         3.5
## 29          5.2         3.4
## 30          4.7         3.2
## 31          4.8         3.1
## 32          5.4         3.4
## 33          5.2         4.1
## 34          5.5         4.2
## 35          4.9         3.1
## 36          5.0         3.2
## 37          5.5         3.5
## 38          4.9         3.6
## 39          4.4         3.0
## 40          5.1         3.4
## 41          5.0         3.5
## 42          4.5         2.3
## 43          4.4         3.2
## 44          5.0         3.5
## 45          5.1         3.8
## 46          4.8         3.0
## 47          5.1         3.8
## 48          4.6         3.2
## 49          5.3         3.7
## 50          5.0         3.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## sROC 0.1-2 loaded
##              Test          Statistic           p value Result
## 1 Mardia Skewness  0.759503524380438 0.943793240544741    YES
## 2 Mardia Kurtosis 0.0934600553610254 0.925538081956867    YES
## 3             MVN               <NA>              <NA>    YES

Pour considérer la distribution comme multinormale, la p-value de la distortion (Mardia Skewness) et la statistique de Kurtosis (Mardia Kurtosis) doit être égale ou plus élevée que 0.05 (Kormaz, 2019, fiche d’aide de la fonction mvn de R). C’est bien le cas pour les données du tableau setosa_sepal.

Retirons de la matrice de covariance les valeurs et vecteurs propres avec la fonction eigen.

Le premier vecteur propre correspond à la première colonne, et le second à la deuxième. Les coordonnées x et y sont les premières et deuxièmes lignes. Les vecteurs propres ont une longueur unitaire (norme de 1). Ils peuvent être mis à l’échelles à la racine carrée des valeurs propres.

Pour effectuer une translation des vecteurs propres au centre du nuage de point, nous avons besoin du centroïde.

On peut observer que, comme je l’ai mentionné plus haut, les vecteurs propres indiquent l’orientation de la covariance, et les valeurs propres indique la longueur associée à cette direction.

9.4.1.1.2 Biplot

Imaginez un nuage de points en 3D, axes y compris. Vous tournez votre nuage de points pour trouver la perspective en 2D qui fera en sorte que vos données soient les plus dispersées possibles. Avec une lampe de poche, vous illuminez votre nuage de points dans l’axe de cette perspective: vous venez d’effectuer une analyse en composantes principales, et l’ombre des points et des axes sur le mur formera votre biplot.

Pour créer un biplot, on juxtapose les descripteurs (variables) en tant que vecteurs propres, représentés par des flèches, et les objets (observations) en tant que scores, représentés par des points. Les résultats d’une ordination peuvent être présentés selon deux types de biplots (Legendre et Legendre, 2012).

Biplot de corrélation permettant de visualiser les corrélations entre des variables météorologiques.

Deux types de projection sont courramment utilisés.

Biplot de distance. Ce type de projection permet de visualiser la position des objets entre eux et par rapport aux descripteurs et d’apprécier la contribution des descripteurs pour créer les composantes principales. Pour créer un biplot de distance, on projette directement les vecteurs propres (\(U\)) en guise de descripteurs. Pour ce qui est des objets, on utilise les scores de l’ACP (\(F\)). De cette manière,

  1. les distances euclidiennes entre les scores sont des approximations des distances euclidiennes dans l’espace multidimentionnel,
  2. la projection d’un objet sur un descripteur perpendiculairement à ce dernier est une approximation de la position de l’objet sur le descripteur et
  3. la projection d’un descripteur sur un axe principal est proportionnelle à sa contribution pour générer l’axe.

Biplot de corrélation. Cette projection permet d’apprécier les corrélations entre les descripteurs. Pour ce faire, les objets et les valeurs propres doivent être transformés. Pour générer les descripteurs, les vecteurs propres (\(U\)) doivent être multipliés par la matrice diagonalisée de la racine carrée des valeurs propres (\(\Lambda\)), c’est-à-dire \(U \Lambda ^{\frac{1}{2}}\). En ce qui a trait aux objets, on multiplie les scores par (\(F\)) par la racine carrée négative des valeurs propres diagonalisées, c’est-à-dire \(F \Lambda ^{- \frac{1}{2}}\). De cette manière,

  1. tout comme c’est le cas pour le biplot de distance, la projection d’un objet sur un descripteur perpendiculairement à ce dernier est une approximation de la position de l’objet sur le descripteur,
  2. la projection d’un descripteur sur un axe principal est proportionnelle à son écart-type et
  3. les angles entre les descripteurs sont proportionnelles à leur corrélation (et non pas leur proximité).

En d’autres mots, le bilot de distances devrait être utilisé pour apprécier la distance entre les objets et le biplot de corrélation devrait être utilisé pour apprécier les corrélations entre les descripteurs. Mais dans tous les cas, le type de biplot utilisé doit être indiqué.

Le triplot est une forme apparentée au biplot, auquel on ajoute des variables prédictives. Le triplot est utile pour représenter les résultats des ordinations contraignantes comme les analyses de redondance et les analyse de correspondance canoniques.

9.4.1.1.3 Application

Bien que l’ACP puisse être effectuée grâce à des modules de base de R, nous utiliserons le module vegan. Le tableau varechem comprend des données issues d’analyse de sols identifiés par leur composition chimique, leur pH, leur profondeur totale et la profondeur de l’humus publiées dans Väre et al. (1995) et exportées du module vegan.

##      N    P     K    Ca   Mg    S    Al    Fe    Mn   Zn  Mo Baresoil Humdepth  pH
## 1 26.6 36.7 171.4 738.6 94.9 33.8  20.7   2.5  77.6  7.4 0.3     23.0      2.8 2.8
## 2 26.2 61.9 202.2 741.2 86.3 48.6 124.3  23.6  94.5 10.2 0.6     56.9      2.5 2.9
## 3 22.8 50.6 151.7 648.0 64.8 30.2  12.1   2.3 122.9  8.1 0.2     23.7      2.6 2.9
## 4 18.0 64.9 224.5 517.6 59.7 52.9 435.1 101.2  38.0  9.5 1.1     21.3      1.8 2.9
## 5 19.1 26.4  61.1 259.1 37.0 21.4 155.1  81.4  20.6  4.0 0.6      5.8      1.9 3.0

Comme nous l’avons vu précdemment, les données de concentration sont de type compositionnelles. Les données compositionnelles du tableau varechem mériteraient d’être transformées (Aitchison et Greenacre, 2002). Utilisons les log-ratios centrés (clr).

## # A tibble: 5 x 15
##        N      P     K    Ca      Mg      S    Al      Fe     Mn    Zn    Mo    Fv Baresoil
##    <dbl>  <dbl> <dbl> <dbl>   <dbl>  <dbl> <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>
## 1 -1.34  -0.760 0.673  2.07  0.549  -0.762 0.181 -1.60   -0.696 -2.26 -5.39  9.34     40.5
## 2 -1.37  -0.655 0.697  2.01 -0.0269 -0.705 0.547 -1.89   -0.869 -2.35 -5.06  9.67     11.2
## 3 -0.773 -0.988 0.175  1.05 -0.174  -0.960 1.50   0.628  -1.88  -3.09 -5.11  9.62     18.6
## 4 -0.871 -0.623 0.872  1.92 -0.244  -0.757 0.250 -1.77   -0.862 -2.41 -5.26  9.76     29.8
## 5 -1.68  -0.469 0.833  2.57  0.347  -0.910 1.04   0.0773 -1.27  -2.39 -7.49  9.32      7.6
## # ... with 2 more variables: Humdepth <dbl>, pH <dbl>

Effectuons l’ACP. Pour cet exemple, nous standardiserons les données étant données que les colonnes Baresoil, Humedepth et pH ne sont pas à la même échelle que les colonnes des clr.

L’objet vareclr_pca contient l’information nécessaire pour mener notre ACP.

## 
## Call:
## rda(X = vareclr_sc) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total              15          1
## Unconstrained      15          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          PC1    PC2    PC3     PC4     PC5     PC6     PC7     PC8     PC9
## Eigenvalue            7.1523 2.4763 2.1122 0.93015 0.57977 0.48786 0.36646 0.29432 0.19686
## Proportion Explained  0.4768 0.1651 0.1408 0.06201 0.03865 0.03252 0.02443 0.01962 0.01312
## Cumulative Proportion 0.4768 0.6419 0.7827 0.84473 0.88338 0.91590 0.94034 0.95996 0.97308
##                          PC10     PC11     PC12     PC13      PC14
## Eigenvalue            0.15434 0.107357 0.095635 0.042245 0.0042200
## Proportion Explained  0.01029 0.007157 0.006376 0.002816 0.0002813
## Cumulative Proportion 0.98337 0.990527 0.996902 0.999719 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  4.309777 
## 
## 
## Species scores
## 
##              PC1     PC2     PC3      PC4        PC5       PC6
## N         0.1437  0.7606 -0.6792  0.19837  0.1128526 -0.050149
## P         0.8670 -0.3214 -0.2950 -0.22940  0.1437960 -0.042884
## K         0.9122 -0.3857  0.2357  0.03469  0.2737020  0.075717
## Ca        0.9649 -0.3362 -0.2147  0.17757 -0.2188717  0.008051
## Mg        0.8263 -0.2723  0.1035  0.52135 -0.1495399 -0.342214
## S         0.8825 -0.3169  0.3539 -0.21216  0.1176279 -0.191386
## Al       -1.0105 -0.2442  0.2146  0.02674 -0.1005560 -0.043569
## Fe       -1.0338 -0.2464  0.1492  0.13162  0.1512218  0.081571
## Mn        0.9556  0.1041 -0.1256 -0.21300  0.2565831  0.275275
## Zn        0.7763 -0.1031 -0.3123 -0.36493 -0.5665691  0.153089
## Mo       -0.2152  0.8717  0.4065 -0.33643 -0.2134335 -0.167725
## Fv        0.2360  0.5776 -0.8112  0.12736  0.1280097 -0.109737
## Baresoil  0.5147  0.4210  0.4472  0.54980 -0.1438570  0.463148
## Humdepth  0.7455  0.4379  0.5194  0.16493  0.0004757 -0.273056
## pH       -0.5754 -0.5864 -0.5957  0.23408 -0.1517661 -0.056641
## 
## 
## Site scores (weighted sums of species scores)
## 
##            PC1       PC2      PC3      PC4      PC5      PC6
## sit1   0.16862  0.423777  0.46731  0.91175  1.10380  1.06421
## sit2  -0.09705 -0.097482  0.61143 -0.29049  1.14916  0.40622
## sit3   0.02831 -0.795737  0.74176 -0.19097 -2.43337 -0.81762
## sit4   1.39081 -0.354376 -0.19377 -0.45160  0.46020 -0.31446
## sit5   1.30346  0.357866  0.29887  0.76856  0.20913 -0.64145
## sit6   0.43636  0.495037  1.21722  1.18128 -0.98242 -0.74474
## sit7   1.07306  0.467575 -0.32245  0.03717  0.13956 -0.64972
## sit8   0.02545  0.659714 -0.28861 -0.01424  0.47105  0.45173
## sit9   1.42005  0.007356 -0.29000 -0.78474  0.97592 -0.80263
## sit10 -0.50638 -0.220909  1.52981  0.26289  0.42135  0.94054
## sit11  0.45392  0.649297  0.44573 -0.26620 -0.74522 -0.53228
## sit12  0.18623  0.259640  0.89112  0.21096 -0.51393  2.24361
## sit13  1.26264  0.225744 -0.96668 -0.69334  0.61990  0.43312
## sit14 -1.48685  0.739545 -0.20926  1.09256  0.61856 -0.87999
## sit15 -0.50622  1.108685 -2.61287 -1.00433 -1.35383  1.21964
## sit16 -1.28653  0.898663 -0.38778 -0.47556 -0.02449 -0.29419
## sit17 -1.72773  0.476962 -0.48878  0.71156  1.06398 -1.33473
## sit18 -0.82844 -0.296515  1.20315 -1.49821 -0.18330  1.05231
## sit19 -1.00247 -0.609253  0.25185 -0.85420  0.71031  0.14854
## sit20 -0.43405 -0.338912  0.55348 -1.35776 -0.81986 -1.02468
## sit21 -0.05083  0.122645 -0.04611 -0.56047 -0.26151 -0.98053
## sit22  0.17891 -2.315489 -0.69084 -0.19547  0.80628  0.04291
## sit23 -0.46443 -2.592018 -1.21615  1.56359 -0.62334  0.28748
## sit24  0.46316  0.728185 -0.49843  1.89726 -0.80791  0.72671

La deuxième ligne de Importance of components, Proportion Explained, indique la proportion de la variance totale captée successivement par les axes principaux. Le premier axe principal comporte 47.68% de la variance. Le deuxième axe principal ajoutant une proportion de 16,51%, une représentation en deux axes principaux présentent 64.19 % de la variance.

##          PC1          PC2          PC3          PC4          PC5          PC6          PC7 
## 0.4768180610 0.1650859388 0.1408156459 0.0620101490 0.0386511040 0.0325238535 0.0244303815 
##          PC8          PC9         PC10         PC11         PC12         PC13         PC14 
## 0.0196215021 0.0131238464 0.0102890284 0.0071571089 0.0063756951 0.0028163495 0.0002813359

La décision du nombre d’axes principaux à retenir est arbitraire. Elle peut dépendre d’un nombre maximal de paramètre à retenir pour éviter de surdimensionner un modèle (curse of dimensionality, section 11) ou d’un seuil de pourcentage de variance minimal à retenir, par exemple 75%. Ou bien, vous retiendrez deux composantes principales si vous désirez présenter un seul biplot.

L’approche de Kaiser-Guttmann (Borcard et al., 2011) consiste à sélectionner les composantes principales dont la valeur propre est supérieure à leur moyenne.

L’approche du broken stick consiste à couper un bâton d’une longueur de 1 en n tranches. La première tranche est de longueur \(\frac{1}{n}\). La tranche suivante est d’une longueur de la tranche précédente à laquelle on aditionne \(\frac{1}{longueur~restante}\). Puis on place les longueurs en ordre décroissant. On retient les composantes principales dont les valeurs propres cumulées sont plus grandes que le broken stick.

Le graphique du broken stick:

Les approches Kaiser-Guttmann et broken stick suggèrent que les trois premières composantes sont suffisantes pour décrire la dispersion des données.

Examinons les loadings (vecteurs propres) plus en particulier. Dans le langage du module vegan, les vecteurs propres sont les espèces (species) et les scores sont les sites.

##                 PC1        PC2        PC3        PC4           PC5          PC6         PC7
## N         0.1437343  0.7606006 -0.6792046  0.1983670  0.1128526122 -0.050148980 -0.09111164
## P         0.8669892 -0.3213683 -0.2949864 -0.2294036  0.1437959857 -0.042883754  0.26894062
## K         0.9122089 -0.3857245  0.2356904  0.0346904  0.2737019601  0.075717162 -0.21662612
## Ca        0.9648855 -0.3361651 -0.2147486  0.1775746 -0.2188716732  0.008050762  0.03630015
## Mg        0.8263327 -0.2723055  0.1035276  0.5213484 -0.1495399242 -0.342213793  0.04617838
## S         0.8824519 -0.3169039  0.3538854 -0.2121562  0.1176278503 -0.191386377 -0.26825994
## Al       -1.0105173 -0.2441785  0.2145614  0.0267422 -0.1005559874 -0.043569364 -0.22737412
## Fe       -1.0337676 -0.2463987  0.1491865  0.1316173  0.1512218115  0.081571443  0.10553041
## Mn        0.9555632  0.1041030 -0.1256178 -0.2130047  0.2565830557  0.275275174  0.20224538
## Zn        0.7763480 -0.1030878 -0.3122919 -0.3649341 -0.5665691228  0.153089144 -0.12332232
## Mo       -0.2152399  0.8717229  0.4064967 -0.3364279 -0.2134335302 -0.167725160  0.13788948
## Fv        0.2360040  0.5775863 -0.8111953  0.1273582  0.1280096553 -0.109737235 -0.20911147
## Baresoil  0.5147445  0.4209983  0.4472351  0.5497950 -0.1438569673  0.463148072 -0.02103009
## Humdepth  0.7455213  0.4379436  0.5193895  0.1649306  0.0004756685 -0.273056212  0.17061078
## pH       -0.5753858 -0.5863743 -0.5957495  0.2340826 -0.1517660977 -0.056640816  0.19890884
##                  PC8          PC9        PC10         PC11         PC12        PC13
## N        -0.06122008  0.315645453  0.08090232 -0.019251478  0.045420621  0.05020956
## P         0.34111276  0.021124287  0.08756299 -0.045741546  0.145128883 -0.03337551
## K        -0.01641260  0.143099440 -0.08737113  0.183005607  0.002260341 -0.10566808
## Ca        0.04775616 -0.073609828 -0.10601799  0.161460554  0.041210064  0.14341793
## Mg       -0.12098602 -0.051599273  0.18373857 -0.009862571 -0.063493608 -0.03782662
## S         0.15822845  0.038378858  0.05100717 -0.138785063 -0.117144869  0.06075094
## Al        0.10598673  0.040586196 -0.14473132 -0.089462074  0.058212507  0.01983102
## Fe       -0.09254655 -0.079426433  0.09908706 -0.006376211  0.049837173 -0.01169516
## Mn       -0.19347804 -0.038859808 -0.07637994 -0.083300112 -0.133353213  0.02679781
## Zn       -0.14862229  0.024026151  0.02643462 -0.064973307  0.051057277 -0.06538348
## Mo        0.17165900  0.032981311  0.01419924  0.128814989 -0.114803631 -0.01989539
## Fv        0.11289753 -0.281443886 -0.08391004 -0.012456867 -0.020157331 -0.05448619
## Baresoil  0.23028292  0.004554036  0.02604286 -0.061147847 -0.019696758 -0.01640490
## Humdepth -0.11310394  0.027515405 -0.23068827 -0.102189307  0.109293684 -0.02485030
## pH        0.12152266  0.150118818 -0.15240317 -0.037691048 -0.153813168 -0.04523353
##                  PC14
## N         0.002340519
## P        -0.010109130
## K         0.001169065
## Ca        0.007419161
## Mg       -0.023575986
## S         0.025874035
## Al       -0.037901576
## Fe        0.036048221
## Mn       -0.021373612
## Zn        0.010896560
## Mo       -0.001335923
## Fv        0.005707928
## Baresoil  0.003823725
## Humdepth  0.016559206
## pH        0.014193061
## attr(,"const")
## [1] 4.309777

L’ordre d’importance des vecteurs propres est établi en ordre croissant des élément des vecteurs propres associées. Un vecteur propre est une combinaison linéaire des variables. Par exemple, le premier vecteur propre pointe surtout dans la direction du Fe (-1.497) et de l’Al (-1.463). Le deuxième pointe surtout vers le Mo (2.145). Les vecteurs (loadings) d’un biplot de distance présentant les des deux premières composantes principales prendront les coordonnées des deux premières colonnes. Le vecteur Al aura la coordonnée [-1.463 ; -0.601], le vecteur de Fe sera placé à [-1.497 ; -0.606] et le vecteur Mo à [-0.312 ; 2.145]. Il existe différentes fonctions d’affichage des biplots. Notez que leur longueur peut être magnifiée pour améliorer la visualisation.

Lançons la fonction biplot pour créer un biplot de distance et un autre de corrélation.

Le biplot de distance permet de dégager les variables qui expliquent davantage la variabilité dans notre tableau: les clr du Fe et de l’Al forment en grande partie le premier axe principal, alors que le clr du Mo forme en grande partie le second axe. Le biplot de corrélation montre que les clr du Fe et du Al sont corrélés dans le même sens, mais das le sens contraire du clr du Mn. L’information sur la teneur en Fe et celle de l’Al est en grande partie redondante. Toutefois, le clr du Mo est presque indépendant du clr du Fe, ceux-ci étant à angle presque droit (~90°). Ces relations peuvent être explorées directement.

Nous avons mentionné que l’ACP est une rotation. Prenons un second exemple pour bien en saisir les tenants et aboutissants. Le tableau de données que nous chargerons provient d’un infographie d’un dauphin, intitullée Bottlenose Dolphin, conçu par l’artiste Tarnyloo. Les points correspondent à la surface d’un dauphin. J’ai ajouté une colonne anatomy, qui indique à quelle partie anatomique le point appartient.

## Parsed with column specification:
## cols(
##   x = col_double(),
##   y = col_double(),
##   z = col_double(),
##   anatomy = col_character()
## )
## # A tibble: 5 x 4
##         x      y        z anatomy   
##     <dbl>  <dbl>    <dbl> <chr>     
## 1 -0.440   1.06   1.24    Caudal fin
## 2  0.0360 -0.494 -0.198   Head      
## 3 -0.177  -0.324  0.00479 Head      
## 4  0.0921 -0.436 -0.0795  Head      
## 5 -0.333   1.07   1.23    Caudal fin

Voici en vue isométrique ce en quoi consiste ce nuage de points.

Effectuons l’ACP sur le dauphin.

On n’y voit pas grand chose, mais si l’on extrait les scores et que l’on raccourcit les vecteurs:

##           PC1         PC2
## x -0.02990131  0.01608095
## y -7.13731672 -1.43221776
## z -4.56612084  2.23859843
## attr(,"const")
## [1] 9.089026

La meilleure représentation du dauphin en 2D, selon la variance, est son profil - en effet, il est plus long et haut que large.

Note. Une ACP effectue seulement une rotation des points. Les distances euclidiennes entre les points sont maintenues.
Note. L’ACP a été conçue pour projetter en un nombre moindre de dimensions des observations dont les distributions sont multinormales (ce n’est évidemment pas le cas du dauphin).
Note. Les axes principaux d’une ACP sont des variables aléatoires. Elles peuvent être assujetties à des tests ststistiques, des modèles, du partitionnement de données, etc.
Excercice. Effectuez maintenant une ACP avec les données d’iris.

9.4.1.2 Analyse de correspondance (AC)

L’analyse de correspondance (AC) est particulièrement appropriée pour traiter des données d’abondance et d’occurence. Tout comme l’analyse en composantes principales, les données apportés vers une AC doivent être dimensionnellement homogènes, c’est-à-dire que chaque variable doit être de même métrique: pour des données d’abondance, cela signifie que les décomptes réfèrent tous au même concept: individus, colonies, surfaces occupées, etc. Alors que la distance euclidienne est préservée avec l’ACP, l’AC préserve la distance du \(\chi^2\), qui est insensible aux double-zéros.

L’AC produit \(min(n,p)-1\) axes principaux orthogonaux qui captent non pas le maximum de variance, mais la proportion de mesures aux carré par rapport à la somme des carrés de la matrice. Le biplot obtenu peut être présenté sous forme de biplot de site (scaling 1), où la distance du \(\chi^2\) est préservée entre les sites ou biplot d’espèces (scaling 2), ou la distance du \(\chi^2\) est préservée entre les espèces. L’AC hérite du coup une propriété importate de la distance du \(\chi^2\), qui accorde davantage de distance entre un compte de 0 et de 1 qu’entre 1 et 2, et davantage entre 1 et 2 qu’entre 2 et 3.

Par exemple, sur ces trois sites, on a compté un individu A de moins que d’individu B.

## # A tibble: 3 x 3
##   Site       A     B
##   <chr>  <dbl> <dbl>
## 1 Site 1     0     1
## 2 Site 2     1     2
## 3 Site 3     9    10

Pourtant, la distance du \(\chi^2\) est plus élevée entre le site 1 et le site 2 qu’entre le site 2 et le site 3.

##           1         2
## 2 0.6724111          
## 3 0.9555316 0.2831205

La distance du \(\chi^2\) donne davantage d’importance aux espèces rares, ce dont une analyse doit tenir compte. Il pourrait être envisageable de retirer d’un tableau des espèces rare, ou bien prétransformer des données d’abondance par une transformation de chord ou de Hellinger (tel que discuté au chapitre 6), puis procéder à une ACP sur ces données (Legendre et Gallagher, 2001).

9.4.1.2.1 Application

Le tableau varespec comprend des données de surface de couverture de 44 espèces de plantes en lien avec les données environnementales du tableau varechem. Ces données ont été publiées dans Väre et al. (1995) et exportées du module vegan.

##   Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube Vacculig Diphcomp
## 1     0.03     3.65     0.00     0.00     4.43     0.00     0.00        0     1.65     0.50
## 2     0.00     1.63     0.35    18.27     7.13     0.05     0.40        0     0.20     0.00
## 3     3.75     5.65     0.00     0.08     5.30     0.10     0.00        0     0.00     0.00
## 4     3.40     0.63     0.00     0.00     1.98     0.05     0.05        0     0.03     0.00
## 5     0.00     5.30     0.00     0.00     8.20     0.00     0.05        0     8.10     0.28
##   Dicrsp Dicrfusc Dicrpoly Hylosple Pleuschr Polypili Polyjuni Polycomm Pohlnuta Ptilcili
## 1    0.0     0.55     0.00     0.00     0.05        0     0.00     0.00     0.03     0.03
## 2    0.3     0.52     0.20     9.97    70.03        0     0.08     0.00     0.07     0.03
## 3    0.0    11.32     0.00     0.00     7.75        0     0.30     0.02     0.07     0.00
## 4    0.0     0.20     0.00     0.00     1.53        0     0.10     0.00     0.05     0.00
## 5    0.0     0.45     0.03     0.00     0.10        0     0.25     0.00     0.03     0.00
##   Barbhatc Cladarbu Cladrang Cladstel Cladunci Cladcocc Cladcorn Cladgrac Cladfimb Cladcris
## 1        0     8.80    29.50    55.60     0.25     0.08     0.25     0.25     0.15     0.10
## 2        0     0.17     0.87     0.00     0.05     0.02     0.03     0.07     0.10     0.02
## 3        0    17.45     1.32     0.12    23.68     0.22     0.50     0.15     0.23     0.97
## 4        0    15.73    20.03    28.20     0.73     0.10     0.15     0.13     0.10     0.15
## 5        0    35.00    42.50     0.28     0.35     0.08     0.20     0.25     0.18     0.13
##   Cladchlo Cladbotr Cladamau Cladsp Cetreric Cetrisla Flavniva Nepharct Stersp Peltapht
## 1     0.03     0.00     0.03   0.00     0.05     0.00     0.15     0.15   0.28        0
## 2     0.00     0.02     0.00   0.00     0.00     0.02     0.00     0.00   0.02        0
## 3     0.00     0.00     0.00   0.00     0.68     0.02     0.00     0.00   0.33        0
## 4     0.00     0.00     0.00   0.05     0.28     0.05    10.03     0.00   0.95        0
## 5     0.08     0.00     0.00   0.00     0.05     0.00     0.23     0.20   0.93        0
##   Icmaeric Cladcerv Claddefo Cladphyl
## 1     0.00     0.00     0.08     0.00
## 2     0.00     0.00     0.08     0.00
## 3     0.02     0.00     1.57     0.05
## 4     0.00     0.05     0.08     0.00
## 5     0.03     0.00     0.10     0.00

Pour effectuer l’AC, nous utiliserons, comme pour l’ACP, le module vegan mais cette fois-ci avec la fonction cca. L’AC en scaling 1 est effectuée sur le tableau des abondances avec les espèces comme colonnes et les sites comme lignes. Les matrices d’abondance transposées indique les sites où chque espèce ont été dénombrées: pour une analyse en scaling 2, on effectue une analyse de correspondance sur la matrice d’abondance (ou d’occurence) transposée.

Pour chacune des AC, je filtre pour m’assurer que toutes les lignes contiennent au moins une observation. Ce n’est pas nécessaire dans notre cas, mais je le laisse pour l’exemple.

## 
## Call:
## cca(X = varespec %>% filter(rowSums(.) > 0)) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.083          1
## Unconstrained   2.083          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7     CA8     CA9
## Eigenvalue            0.5249 0.3568 0.2344 0.19546 0.17762 0.12156 0.11549 0.08894 0.07318
## Proportion Explained  0.2520 0.1713 0.1125 0.09383 0.08526 0.05835 0.05544 0.04269 0.03513
## Cumulative Proportion 0.2520 0.4233 0.5358 0.62962 0.71489 0.77324 0.82868 0.87137 0.90650
##                          CA10    CA11    CA12    CA13     CA14     CA15     CA16     CA17
## Eigenvalue            0.05752 0.04434 0.02546 0.01710 0.014896 0.010160 0.007830 0.006032
## Proportion Explained  0.02761 0.02129 0.01222 0.00821 0.007151 0.004877 0.003759 0.002896
## Cumulative Proportion 0.93411 0.95539 0.96762 0.97583 0.982978 0.987855 0.991614 0.994510
##                           CA18     CA19      CA20      CA21      CA22      CA23
## Eigenvalue            0.004008 0.002865 0.0019275 0.0018074 0.0005864 0.0002434
## Proportion Explained  0.001924 0.001375 0.0009253 0.0008676 0.0002815 0.0001168
## Cumulative Proportion 0.996434 0.997809 0.9987341 0.9996017 0.9998832 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                 CA1       CA2      CA3       CA4        CA5       CA6
## Callvulg  0.0303167 -1.597460  0.11455 -2.894569  0.1376073  2.291129
## Empenigr  0.0751030  0.379305  0.39303  0.023675  0.8568729 -0.400964
## Rhodtome  1.1052309  1.499299  3.04284  0.120106  3.2324306 -0.283510
## Vaccmyrt  1.4614812  1.622935  2.72375  0.231688  0.4604556  0.712538
## Vaccviti  0.1468014  0.313436  0.14696  0.243505  0.6868371 -0.147815
## Pinusylv -0.4820096  0.588517 -0.36020 -0.127094  0.4064754  0.386604
## Descflex  1.5348239  1.218806  1.87562 -0.001340 -1.3136979 -0.070731
## Betupube  0.6694503  1.951826  3.84017  1.389423  7.5959115 -0.244478
## Vacculig -0.0830789 -1.629259  1.05063  0.802648 -0.3058811 -1.625341
## Diphcomp -0.5446464 -1.037570  0.52282  0.940275  0.3682126 -1.082929
## Dicrsp    1.8120408  0.360290 -4.92082  3.088562  1.3867372  0.157815
## Dicrfusc  1.2704743 -0.562978 -0.39718 -2.929542  0.3848272 -2.408710
## Dicrpoly  0.7248118  1.409347  0.80341  1.915549  4.5674148  1.295447
## Hylosple  2.0062408  1.743883  2.27549  0.928884 -3.7648428  2.254851
## Pleuschr  1.3102086  0.583036 -0.01004  0.137298 -1.1216144  0.200422
## Polypili -0.3805097 -1.243904  0.54593  1.477188 -0.7276341 -0.387641
## Polyjuni  1.0133795  0.099043 -2.24697  1.510641  0.7729714 -3.062378
## Polycomm  0.8468241  1.321773  1.13585  1.140723  2.6836594 -0.605038
## Pohlnuta -0.0136453  0.589290 -0.35542  0.135481  0.9369707  0.397246
## Ptilcili  0.4223631  1.598584  3.43474  1.400065  6.3209491  0.198935
## Barbhatc  0.5018348  2.119334  4.57303  1.693188  8.1101807  0.645995
## Cladarbu -0.1531729 -1.483884  0.20024  0.193680  0.0734141  0.358926
## Cladrang -0.5502561 -1.084008  0.40552  0.724060 -0.3357992 -0.335924
## Cladstel -1.4373146  1.077753 -0.44397 -0.375926 -0.2421525  0.004212
## Cladunci  0.8151727 -1.006186 -1.82587 -1.389523  1.6046713  3.675908
## Cladcocc -0.2133215 -0.584429 -0.21434 -0.567886 -0.0003788 -0.145303
## Cladcorn  0.2631227 -0.177858 -0.44464  0.272422  0.3992282 -0.306738
## Cladgrac  0.1956947 -0.311167 -0.23894  0.379013  0.4933026  0.037581
## Cladfimb  0.0009213 -0.161418  0.18463 -0.435908  0.4831233 -0.143751
## Cladcris  0.3373031 -0.470369 -0.05093 -0.823855  0.7182250  0.636140
## Cladchlo -0.6200021  1.207278  0.21889  0.426447  1.9506082  0.120722
## Cladbotr  0.5647242  1.047333  2.65330  0.907734  4.4946805  1.201655
## Cladamau -0.6598144 -1.512880  0.83251  1.577699 -0.0407227 -1.419139
## Cladsp   -0.8209003  0.476164 -0.49752 -0.998241 -0.2393208  0.390785
## Cetreric  0.2458192 -0.689228 -1.68427 -0.131681  0.7439412  2.374535
## Cetrisla -0.3465221  1.362693  0.85897  0.396752  2.7526968  0.396591
## Flavniva -1.4391907 -0.833589 -0.12919  0.007071 -1.4841375  2.956977
## Nepharct  1.6813309  0.199484 -4.33509  2.229917  0.9561223 -5.472858
## Stersp   -0.5172793 -2.280900  0.99775  2.377013 -0.8892757 -1.441228
## Peltapht  0.4035858 -0.043265  0.04538  0.711040  0.1824679 -0.841227
## Icmaeric  0.0378754 -2.419595  0.72135  0.361302 -0.3736424 -2.092136
## Cladcerv -0.9232858 -0.005233 -1.22058  0.305290 -0.8142627  0.414135
## Claddefo  0.5190399 -0.496632 -0.15271 -0.695927  0.9042143  0.909191
## Cladphyl -1.2836161  1.155872 -0.79912 -0.741170 -0.1608002  0.490526
## 
## 
## Site scores (weighted averages of species scores)
## 
##             CA1      CA2       CA3      CA4        CA5      CA6
## sit1  -0.108122 -0.53705  0.229574  0.24412  0.1405624 -0.14253
## sit2   0.697118 -0.14441 -0.031788 -0.21743 -0.2738522 -0.08146
## sit3   0.987603  0.15042 -1.348447  0.80472  0.3095168  0.46773
## sit4   0.851765  0.49901  0.443559  0.12277 -0.4814871  0.07589
## sit5   0.359881 -0.05608  0.145813  0.15087  0.2405263 -0.17770
## sit6   0.003545  0.37017  0.027760  0.06168 -0.1158930 -0.03413
## sit7   0.860732 -0.11504  0.110869 -1.02169  0.0772348 -0.60530
## sit8   0.636936 -0.33250  0.001120 -0.79797  0.0130769 -0.54049
## sit9   1.279352  0.81557  0.670053  0.23137 -0.8929976  0.41783
## sit10 -0.195009 -0.80564  0.117686 -0.58286 -0.0007212  0.53071
## sit11  0.528532 -0.70420 -0.517771 -0.86836  0.5713441  0.91671
## sit12  0.382866 -0.18686 -0.004789  0.10156  0.0458125  0.21087
## sit13  0.990715  0.11967 -1.110040  0.44929  0.1885902 -0.70694
## sit14 -0.264704 -1.06013  0.334900  0.45973 -0.0326631 -0.19945
## sit15 -0.428410 -1.20765  0.374344  0.74970 -0.2596294 -0.30467
## sit16 -0.330534 -0.77498  0.130760  0.22391  0.0632686  0.09060
## sit17 -0.899601  0.12075 -0.075742  0.03842 -0.1489585 -0.12031
## sit18 -0.770294 -0.35351 -0.033779 -0.01795 -0.3007839  0.44303
## sit19 -0.992193  0.50319 -0.157505 -0.07070 -0.1065172 -0.09928
## sit20 -0.937173  0.78688 -0.258119 -0.19377 -0.0343535 -0.01259
## sit21 -0.726413  0.49163 -0.157235 -0.08698 -0.0105774 -0.02801
## sit22 -1.002083  0.71239 -0.236526 -0.18643 -0.0231666 -0.04928
## sit23 -0.322647 -0.03871 -0.001297  0.09029 -0.1481448  0.06934
## sit24  0.259527  0.80746  1.124258  0.36083  1.5437866  0.07051

Créons les biplots.

Le biplot des espèces, à gauche (scaling = 1), montre la distribution des sites selon les espèces. Les emplacements des scores (en noir) montrent les contrastes entre sites selon les espèces qui les recouvrent. Les sites 14 et 15, par exemple, contrastent les sites 19, 20, 21 et 22 selon le 2ième axe principal. Par ailleurs, les axes principaux sont formé de plusieurs espèces dont aucune ne domine clairement.

Le biplot des sites, à droite (scaling = 2), montre la distribution des recouvrements d’espèces selon les sites. Par exemple, les espèces Betupube (Betula pubescens) et Barbhatc (Barbilophozia hatcheri ) se recouvrent en particulier le site 24. Le site 1 est difficile à identifier, car il est couvert par plusieurs noms d’espèces, au bas au centre. Les sites 3 et 13 se confondent avec Dicrsp (une espèce de Dicranum) qui le recouvre amplement.

Pour les deux types de biplot, les sites où les espèces situés près de l’origine, car ils peuvent être soit près de la moyenne, soit distribués uniformément.

Le nombre de composantes à retenir peut être évalué par les approches Kaiser-Guttmann et broken-stick.

Pour les deux scalings, l’approche Kaiser-Guttmann propose 7 axes, tandis que l’approche broken-stick en propose 5.

Les représentations biplot d’analyse de correspondance peuvent prendre la forme d’un boomerang, en particulier celles qui sont basées sur des données d’occurence. Le tableau suivant initialement de Chessel et al. (1987) et est distribué dans le module ade4.

Les numéros de sites correspondent à la position dans une rivière, 1 étant en amont et 30 en aval. Le premier axe discrimine l’amont et l’aval, tandis que le deuxième montre deux niches en amont. Bien que l’on observe une discontinuité dans le cours d’eau, il y a une continuité dans les abondances. Cet effet peut être corrigé en retirant la tendance de l’analyse de correspondance par une detrended correspondance analysis. Pour cela, il faudra utiliser la fonction decorana, ce qui ne sera pas couvert ici.

L’analyse des correspondances multiples (ACM) est utile pour l’ordination des données catégorielles. Le module ade4 est en mesure d’effectuer des AMC, mais n’est pas couvert dans ce manuel.

Excercice. Effectuez et analysez une AC avec les données de recouvrement varespec.

9.4.1.3 Positionnement multidimensionnel (PoMd)

Le positionnement multidimensionnel (PoMd), ou manifold analysis, se base sur les assiciations entre les objets (mode Q) ou les variables (mode R) pour en réduire les dimensions. Alors que l’analyse en composantes principales conserve la distance euclidienne et que l’analyse de correspondance conserve la distance du \(\chi^2\), le PoMd conserve l’association que vous sélectionnerez à votre convenance. Le PoMd vise à représenter en un nombre limité de dimensions (souvent 2) la distance (ou dissimilarité) qu’ont les objets (ou des variables) les uns par rapport aux autres dans l’espace multidimensionnel.

Il existe deux types d’AEM. Le PoMd-métrique (metric multidimentional scaling MMDS, parfois le metric est retiré, MDS, et parfois l’on parle de classic MDS) vise à représenter fidèlement la distance entre les objets ou les variables. Le PoMd-métrique ne devrait être utilisée que lorsque la métrique n’est ni euclidienne, ni de \(\chi^2\) et que l’on désire préserver les distances entre les objets. L’PoMd-métrique aussi appelée analyse en coordonnées principales (ACoP ou de l’anglais PCoA) .

Le PoMd-non-métrique (nonmetric multidimentional scaling, NMDS) vise quant à lui à représenter l’ordre des distances entre les objets ou les variables. C’est une approche par rang: le PoMd-non-métrique vise représenter les objets sont plus proches ou plus éloignées les uns des autres plutôt que de représenter leur similarité dans l’espace multidimentionnelle.

L’IsoMap, pour isometric feature mapping, est une extension du PoMd qui recontruit les distances selon les points retrouvés dans le voisinage. Les isomaps sont en mesure d’applatir des données ayant des formes complexes.

Nous ne traitons pour l’instant que de l’PoMd-métrique (fonction vegan::cmdscale) et des PoMd-non-métrique (fonction vegan::metaMDS).

9.4.1.3.1 Application

Utilisons les données d’abondance que nous avions au tout début de ce chapitre. La matrice d’association de Bray-Curtis sera utilisée.

Les sites 2 et 3 devraient être plus près l’un et l’autre, puis les sites 3 et 4. Les autres associations sont éloignés d’environ la même distance. Lançons le calcul de la PoMd-métrique.

## species scores not available

On observe en effet que les sites 2 et 3 sont les plus près. Les sites 3 et 4sont plus éloignés. Les sites 1, 2 et 4 font à peu près un triangle équilatéral, ce qui correspond à ce à quoi on devrait s’attendre. Les wa-scores permettent de juxtaposer les espèces sur les sites, pour référence. Le colibri n’est présent que sur le site 2. Le site 1 est populé par des jaseurs et des mésanges, et c’est le seul site où l’on a observé une citelle. On a observé des chardonnerets sur les sites 2 et 3. Sur le site 4, on n’a observé que des bruants, que l’on a aussi observé ailleurs, sauf au site 2.

Le PoMd-non-métrique (non metric dimensional scaling, NMDS) fonctionne de la même manière que la PoMd-métrique, à la différence que la distance est basée sur les rangs. À cet égard, le site 4 à une distance de 0.76 du site 3, mais plutôt le deuxième plus loin, après le site 2 et avant le site 1. Utilisons la fonction metaMDS.

## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.1712075  max resid 0.2311738 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1038622  max resid 0.1447403 
## Run 3 stress 0 
## ... Procrustes: rmse 0.1355982  max resid 0.1807723 
## Run 4 stress 0 
## ... Procrustes: rmse 0.1149604  max resid 0.1456473 
## Run 5 stress 0 
## ... Procrustes: rmse 0.08167542  max resid 0.1146596 
## Run 6 stress 0 
## ... Procrustes: rmse 0.1284115  max resid 0.1569041 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1834823  max resid 0.2366359 
## Run 8 stress 0 
## ... Procrustes: rmse 0.1367236  max resid 0.1848335 
## Run 9 stress 0 
## ... Procrustes: rmse 0.09492646  max resid 0.1233517 
## Run 10 stress 0 
## ... Procrustes: rmse 0.12845  max resid 0.1738103 
## Run 11 stress 8.215235e-05 
## ... Procrustes: rmse 0.07684147  max resid 0.1148002 
## Run 12 stress 0 
## ... Procrustes: rmse 0.06054277  max resid 0.07388255 
## Run 13 stress 0 
## ... Procrustes: rmse 0.1239983  max resid 0.1544609 
## Run 14 stress 0 
## ... Procrustes: rmse 0.09964013  max resid 0.1300011 
## Run 15 stress 0 
## ... Procrustes: rmse 0.06571301  max resid 0.08654521 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1092868  max resid 0.1707372 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1520543  max resid 0.202858 
## Run 18 stress 0 
## ... Procrustes: rmse 0.158808  max resid 0.2101473 
## Run 19 stress 0 
## ... Procrustes: rmse 0.1124983  max resid 0.1450684 
## Run 20 stress 0 
## ... Procrustes: rmse 0.09344561  max resid 0.1301929 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(assoc_mat, k = nrow(abundance) - 1, eig = TRUE): stress is (nearly) zero:
## you may have insufficient data
## species scores not available

Dans ce cas, entre PoMd-métrique et non-métrique, les résultats peuvent être interprétés de manière similaire.

En ce qui a trait au dauphin,

Pour plus de détails, je vous invite à vous référer à Borcard et al. (2011)) ou de consulter l’excellent site GUSTA ME.

9.4.1.4 Conclusion sur l’ordination non contraignante

Lorsque les données sont euclidiennes, l’analyse en composantes principales (ACP) dervait être utilisée. Lorsque la métrique est celle du \(\chi^2\), on préférera l’analyse de correspondance (AC). Si la métrique est autre, le positionnement multidimensionel (PoMd) est préférable. Dans ce dernier cas, si l’on recherche une représentation simplifiée de la distance entre les objets ou variables, on utilisera un PoMd-métrique. À l’inverse, si l’on désire une représentation plus fidèle au rang des distances, on préférera l’PoMd-non-métrique.

9.4.2 Ordination contraignante

Alors que l’ordination non contraignante vous permet de dresser un protrait de vos variables, l’ordination contraignante (ou canonique) permet de tester statistiquement ainsi que de représenter la relation entre plusieurs variables explicatives (par exemple, des conditions environnementales) et une ou plusieurs variables réponses (par exemple, les espèces observées).

  • L’analyse discriminante n’a fondamentalement qu’une seulement variable réponse, et celle-ci doit décrire l’appartenance à une catégorie.
  • L’analyse de redondance sera préférée lorsque le nombre de variable est plus restreint (variables ionomiques et indicateurs de performance des cultures). Les détails, ainsi que les tenants et aboutissants de ces méthodes, sont présentés dans Numerical Ecology (Legendre et Legendre, 2012).
  • L’analyse canonique des corrélations sera préférée lorsque les variables sont parsemées (beaucoup de colonnes avec beaucoup de zéros, comme les variables d’abondance).

9.4.2.1 Analyse discriminante

Alors que l’analyse en composante principale vise à présenter la perspective (les axes) selon laquelle les points sont les plus éclatées, l’analyse discriminante, le plus souvent utilisé dans sa forme linéaire (ADL) et quadratique (ADQ), vise à présenter la perspective selon laquelle les groupes sont les plus éclatés, les groupes formant la variable contraignante. Ces groupes peuvent être connus (e.g. cultivar, région géographique) ou attribués (exemple: par partitionnement). L’ADL est parfois nommée analyse canonique de la variance.

L’AD vise à représenter des différences entre des groupes aux moyens de combinaisons linéaires (ADL) ou quadratique (ADQ) de variables mesurées. Sa représentation sous forme de biplot permet d’apprécier les différences entre les groupes d’identifier les variables qui sont responsables de la discrimination.

Biplot de distance de l’analyse discriminante des ionomes d’espèces de plantes à fruits cultivées sauvages et domestiquées, Source: Parent et al. (2013)

L’ADL a été développée par Fisher (1936), qui à titre d’exemple d’application a utilisé un jeu de données de dimensions d’iris collectées par Edgar Anderson, du Jardin botanique du Missouri, sur 150 spécimens d’iris collectés en Gaspésie (Est du Québec), ma région natale (suis-je assez chauvin?). Ce jeu de données est amplement utilisé à titre d’exemple en analyse multivariée.

Williams (1983) a présenté les tenants et aboutissants de l’ADL en écologie. Tout comme les données passant pas une ACP doivent suivre une distribution multinormale pour être statistiquement valide, les distributions des groupes dans une ADL doivent être multinormales et les variances des points par groupe doivent être homogènes… ce qui est rarement le cas en science. Néanmoins:

Heureusement, il y a des évidences dans la littérature que certaines d’entre [ces règles] peuvent être transgressées modérément sans de grands changement dans les taux de classification. Cette conclusion dépends, toutefois, de la sévérité des transgressions, et de facteurs structueaux comme la position relative des moyennes des populations et de la nature des dispersions. - Williams (1983)

L’ADL peut servir autant d’outil d’interprétation que d’outil de classification, c’est à dire de prédire une catégorie selon les variables (chapitre 11). Dans les deux cas, lorsque le nombre de variables approchent le nombre d’observation, les résultats d’une ADL risque d’être difficilement interprétables. Le test approprié pour évaluer l’homodénéité de la covariance est le M-test de Box. Ce test est peu documenté dans la littérature, est rarement utilisé mais a la réputation d’être particulièrement sévère.

Il est rare que des données écologiques aient des dispersions (covariances) homogènes. Contrairement à l’ADL, l’ADQ ne demande pas à ce que les dispersions (covariances) soient homogènes. Néanmoins, l’ADQ ne génère ni de scores, ni de loadings: il s’agit d’un outil pour prédire des catégories (classification), non pas d’un outil d’ordination.

9.4.2.1.1 Application

Utilisons les données d’iris.

Testons la multinormalité par groupe. Rappelons-nous que pour considérer la distribution comme multinormale, la p-value de la distortion ainsi que la statistique de Kurtosis doivent être égale ou plus élevée que 0.05. La fonction split sépare le tableau en listes et la fonction map applique la fonction spécifiée à chaque élément de la liste. Cela permet d’effectuer des tests de multinormalité sur chacune des espèces d’iris.

## $setosa
##              Test        Statistic           p value Result
## 1 Mardia Skewness 25.6643445196298 0.177185884467652    YES
## 2 Mardia Kurtosis 1.29499223711605 0.195322907441935    YES
## 3             MVN             <NA>              <NA>    YES
## 
## $versicolor
##              Test         Statistic           p value Result
## 1 Mardia Skewness  25.1850115362466 0.194444483140265    YES
## 2 Mardia Kurtosis -0.57186635893429 0.567412516528727    YES
## 3             MVN              <NA>              <NA>    YES
## 
## $virginica
##              Test         Statistic           p value Result
## 1 Mardia Skewness  26.2705981752915 0.157059707690356    YES
## 2 Mardia Kurtosis 0.152614173978342 0.878702546726567    YES
## 3             MVN              <NA>              <NA>    YES

Le test est passé pour toutes les espèces. Voyons maintenant l’homogénéité de la covariance. Pour ce faire, nous aurons besoin de la fonction boxM, disponible avec le module biotools. Pour que les covariances soient considérées comme égales, la p-vaule doit être supérieure à 0.05.

## 
## Attaching package: 'heplots'
## The following object is masked from 'package:pls':
## 
##     coefplot
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  iris %>% select(-Species)
## Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16

On est loin d’un cas où les distributions sont homogènes. Nous allons néanmoins procéder à l’analyse discriminante avec le module ade4. Nous aurons d’abord besoin d’effectuer une ACP avec la fonction dudi.pca de ade4 (en spécifiant une mise à l’échelle), que nous projeterons en ADL avec discrimin.

La visualisation peut être effectuée directement sur l’objet issu de la fonction discrimin.

Il s’agit toutefois d’une visualisation pour le diagnostic davantage que pour la publication. Si l’objectif est la pubilcation, vous pourriez utiliser la fonction plotDA que j’ai conçue à cet effet. J’ai aussi conçu une fonction similaire qui utilise le module graphique de base de R.

## Loading required package: ellipse
## 
## Attaching package: 'ellipse'
## The following object is masked from 'package:car':
## 
##     ellipse
## The following object is masked from 'package:graphics':
## 
##     pairs
## Loading required package: grid
## Loading required package: plyr
## ---------------------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ---------------------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:plotly':
## 
##     arrange, mutate, rename, summarise
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise, summarize
## The following object is masked from 'package:purrr':
## 
##     compact

À la différence de l’ACP, l’ADL maximise la sépatation des groupes. Nous avions noté avec l’ACP que les dimensions des pétales distingaient les groupes. Puisque nous avions justement des informations sur les groupes, nous aurions pu procéder directement à un ADL pour obtenir des conclusions plus directes. Si la longueur des pétales permet de distinguer l’espèce setosa des deux autres, la largeur des pétales permet de distinguer virginica et versicolor, bien que les nuages de points se superposent. De manière bivariée, les régions de confiance des moyennes des scores discriminants (petites ellipses) montrent des différence significatives au seuil 0.05.

Excercice. Si l’on effectuait l’ADL sur notre dauphin, avec la colonne anatomy comme variable de regroupement, qu’obtiendrions-nous? Si l’on consière la nageoire codale (queue) comme faisant partie du corps? Quelles sont les limitations?

9.4.2.2 Analyse de redondance (RDA)

En anglais, on la nomme redundancy analysis, souvent abrégée RDA. Elle est utilisée pour résumer les relations linéaires entre des variables réponse et des variables explicatives. La “redondance” se situe dans l’utilisation de deux tableaux de données contenant de l’information concordante. L’analyse de redondance est une manière élégante d’effectuer une régresssion linéaire multiple, où la matrice de valeurs prédites par la régression est assujettie à une analyse en composantes principales. Il est ainsi possible de superposer les scores des variables explicatives à ceux des variables réponse.

Plus précisément, une RDA effectue les étapes suivantes (Borcard et al. (2011)) entre une matrice de variables indépendantes (explicatives) \(X\) et une matrice de variables dépendantes (réponse) \(Y\).

9.4.2.2.1 1. Régression entre \(Y\) et \(X\)

Pour chacune des variables réponse de \(Y\) (\(y_1\), \(y_2\), , \(y_j\)), effectuer une régression linéaire sur les variables explicatives \(X\).

\[\hat{y}_j = b_j + m_{1, j} \times x_1 + m_{2, j} \times x_2 + ... + m_{i, j} \times x_i\]

\[\hat{y}_j = y_j + y_{res, j}\]

Pour chaque observation (\(n\)), nous obtenons une série de valeurs de \(\hat{y}_j\) et de \(y_{res, j}\). Donc chaque cellule de la matrice \(Y\) a ses pendant \(\hat{y}\) et \(y_{res}\). Nous obtenons ainsi une matrice de prédiction \(\hat{Y}\) et une matrice des résidus \(Y_{res} = Y - \hat{Y}\).

9.4.2.2.2 2. Analyse en composantes principales

Ensuite, on effectue une analyse en composantes principales (ACP) sur la matrice des prédictions \(\hat{Y}\). On obtient ainsi ses valeurs et vecteurs propres. Nommons \(U\) ses vecteurs propres. Les fonctions de RDA mettent souvent ces veceturs à l’échelle avant de les retourner à l’utilisateur. En ordination écologique, ces vecteurs mis à l’échelle sont souvent appelés les scores des espèces, bien qu’il ne s’agisse pas nécessairement d’espèces, mais plus généralement des variables de la matrice dépendante \(Y\).

Il est aussi possible d’effectuer une ACP sur \(Y_{res}\).

9.4.2.2.3 3. Calculer les scores

Les vecteurs propres \(U\) sont utilisés pour calculer les scores des sites, \(Y \times U\), ainsi que les contraintes de site \(\hat{Y} \times U\).

9.4.2.2.4 Application

Nous allons utiliser la fonction rda du module vegan. En ce qui a trait aux données, utilisons les données varespec (matrice Y) et varechem (matrice X). La fonction rda peut fonctionner avec l’interface-formule de R, où à gauche du ~ on retrouve le Y (la matrice de la communauté écologique, i.e. les abondances d’espèces) contre le X (l), à gauche, ce qui peut être pratique pour l’analyse d’intéractions. Mais pour comparer deux matrices, nous pouvons définir X et Y. Ce qui est mélangeant, c’est que vegan, contrairement aux conventions, défini X comme étant la matrice réponse et Y comme étant la matrice explicative.

La fonction ordiplot permet de créer un triplot de base. La représentation des wascores est réputée plus robuste (moins susceptible d’être bruitée), mais leur interprétation porte à confusion (Borcard et al. (2011)).

Triplot de distance (scaling 1). Les angles entre les variables explicatives représentent leur corrélation (non pas les variables réponse).

Triplot de corrélation (scaling 2). Les angles entre les variables représentent leurs corrélation, que les variables soient réponse ou explicative, ou entre variables réponses et variables explicatives. Les distances entre les objets sur le triplot ne sont pas des approximation de leur distance euclidienne.

Les triplots montrent que les variables ont toutes un rôle important sur la dispersion des sites autours des axeds principaux. Le premier axe principal est composé de manière plus marquée par le clr de l’Al et celui du Fe. Le deuxième axe principal est composé de manière plus marquée par le clr du S, du P et du K. Le triplot de corrélation ne présente pas de tendance appréciable pour la plupart des espèces, qui ne possèdent pas de niche particulière. Toutefois, l’espèce Cladstel, présente surtout dans les sites 9 et 10, est liée à de basses teneurs en N et à de faibles valeurs de Baresoil (sol nu). L’espèce Pleuschr est liée à des sols où l’on retrouve une grande épaisseur d’humus, ainsi que des teneurs élevées en nutriment K, P, S, Ca, Mg et Zn. Elle semble apprécier les sols à bas pH, mais à faible teneur en Fe et Al. La teneur en N lui semble plus indifférente (son vecteur étant presque perpendiculaire).

On pourra personnaliser les graphiques en extrayant les scores.

On pourra effectuer une analyse de Kaiser-Guttmann ou de broken-stick de la même manière que précédemment. Étant une collection de régressions, une RDA est en mesure d’effectuer des tests statistiques sur les coefficients de la régression en utilisant des permutations pour tester la signification des coefficients et des axes d’une RDA. On doit néanmoins obligatoirement effectuer la RDA avec l’interface formule. L’a variable de gauche’objet à gauche du ~ peut être une matrice ou un tableau, et celui de droite est défini dans data. Le . dans l’interface formule signifie “une combinaison linéaire de toutes les variables, sans intéraction”.

La signification des axes est difficile à interpréter. Toutefois, celui des variables présente un intérêt.

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Fv + Baresoil + Humdepth + pH, data = vareclr, scale = FALSE)
##          Df Variance      F Pr(>F)   
## N         1   216.13 4.8470  0.009 **
## P         1   272.71 6.1159  0.003 **
## K         1   194.97 4.3724  0.016 * 
## Ca        1    24.92 0.5589  0.686   
## Mg        1    52.61 1.1799  0.303   
## S         1   100.07 2.2441  0.095 . 
## Al        1   177.91 3.9900  0.024 * 
## Fe        1   118.59 2.6595  0.053 . 
## Mn        1    25.96 0.5822  0.621   
## Zn        1    35.81 0.8030  0.483   
## Mo        1    23.51 0.5273  0.668   
## Baresoil  1    98.64 2.2122  0.100 . 
## Humdepth  1    43.59 0.9777  0.402   
## pH        1    38.93 0.8730  0.460   
## Residual  9   401.31                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-value est la probabilité que les pentes calculées pour les variables émergent de distributions dont la moyenne est nulle. Au seuil 0.05, les variables significatives sont (les clr de) l’azote, le phosphore, le potassium et l’aluminium.

Dans le cas des matrices d’abondance (ce n’est pas le cas de varespec, constituée de données de recouvrement), il est préférable avec les RDA de les transformer préalablement avec la transformation compositionnelle, de chord ou de Hellinger (chapitre 8). Une autre option est d’effectuer une RDA sur des matrices d’association en passant par une analyse en coordonnées principales (Legendre et Anderson, 1999). Enfin, les données d’abondance à l’état brutes devraient plutôt passer utiliser une analyse canonique des corrélations.

9.4.2.3 Analyse canonique des correspondances (ACC)

L’analyse canonique des correspondances (Canonical correspondance analysis), ACC, a été à l’origine conçue pour étudier les liens entre des variables environnementales et l’abondance (décompte) ou l’occurence (présence-absence) d’espèces (ter Braak, 1986). L’ACC est à la RDA ce que la CA est à l’ACP. Alors que la RDA préserve les distance euclidiennes entre variables dépendantes et indpendantes, l’ACC préserve les distances du \(\chi^2\). Tout comme l’AC, elle hérite du coup une propriété importate de la distance du \(\chi^2\): il y a davantage davantage d’importance aux espèces rares.

L’analyse des correspondances canoniques est souvent utilisée dans la littérature, mais dans bien des cas une RDA sur des données d’abondance transformées donnera des résultats davantage intérprétables (Legendre et Gallagher, 2001).

9.4.2.3.1 Application

Cet exemple d’application concerne des données d’abondance. Nous allons conséquemment utiliser une CCA avec la fonction cca, toujours avec le module vegan.

Les tableaux doubs_fish et doubs_env comprennent respectivement des données d’abondance d’espèces de poissons et dans différents environnements de la rivière Doubs (Europe) publiées dans Verneaux. (1973) et exportées du module ade4.

Sur le site no 8, aucun poisson n’a pas été observé. Les observations ne comprenant que des zéro doivent être préalablement retirées.

De la même manière qu’avec la fonction rda de vegan, nous utilisons cca pour l’ACC.

Comparons les résultats

Triplot de distance (scaling 1).

  1. La projection des variables réponse à angle droit sur les variables explicatives est une approximation de la réponse sur l’explication. (2) Un objet (site ou réponse) situé près d’une variable explicative est plus susceptible d’avoir le décompte 1. (3) Les distances entre les variables (réponse et explicatives) approximent la distance du \(\chi^2\) (traduction adaptée de Borcard et al. (2011)).

Triplot de corrélation (scaling 2).

  1. La valeur optmiale de l’espèce sur une variable environnementale quantitative peut être obtenue en projetant l’espèce à angle droit sur la variable. (2) Une espèce se trouvant près d’une variable environnementale est susceptible de se trouver en plus grande abondance aux sites de statut 1 pour cette variable. (3) Les distances n’approximent pas la distance du \(\chi^2\) (traduction adaptée de Borcard et al. (2011)).