1 L’analyse compositionnelle

Non-non, Bugs. Nous n’allons pas composer de la musique, mais apprendre à transformer des données qui font partie d’un tout.

Oh non, Bugs. Nous n’allons pas travailler sur la portée, mais bien en langage R. Par ailleurs, les fichiers de calculs qui génèrent ces notes sont disponibles sur GitHub.

1.1 Contexte

En 1898, le statisticien Karl Pearson nota que des corrélations étaient induites lorsque l’on effectuait des ratios par rapport à une variable commune.

Source [Karl Pearson, 1897. Mathematical contributions to the theory of evolution.—on a form of spurious correlation which may arise when indices are used in the measurement of organs. Proceedings of the royal society of London](https://royalsocietypublishing.org/doi/pdf/10.1098/rspl.1896.0076)

Figure 1.1: Source Karl Pearson, 1897. Mathematical contributions to the theory of evolution.—on a form of spurious correlation which may arise when indices are used in the measurement of organs. Proceedings of the royal society of London

Faisons l’exercice! Nous générons au hasard 1000 données (comme le proposait Pearson) pour trois dimensions: le fémur, le tibia et l’humérus. Ces dimensions ne sont pas générées par des distributions corrélées.

library("tidyverse")
set.seed(3570536)
n <- 1000
bones <- tibble(femur = rnorm(n, 10, 3),
                tibia = rnorm(n, 8, 2),
                humerus = rnorm(n, 6, 2))
plot(bones)
Représentation bivariées de longueur fictives des os, générées au hasard.

Figure 1.2: Représentation bivariées de longueur fictives des os, générées au hasard. 

cor(bones)
##                femur        tibia      humerus
## femur    1.000000000 -0.069006171  0.002652292
## tibia   -0.069006171  1.000000000 -0.008994704
## humerus  0.002652292 -0.008994704  1.000000000

Pourtant, si j’utilise des ratios allométriques avec l’humérus comme base,

bones_r <- bones %>% 
  transmute(fh = femur/humerus,
            th = tibia/humerus)
plot(bones_r)
text(30, 20, paste("corrélation =", round(cor(bones_r$fh, bones_r$th), 2)), col = "blue")
Corrélation entre relations allométriques d'attributs générés au hasard.

Figure 1.3: Corrélation entre relations allométriques d’attributs générés au hasard.

Pas de panique, Bugs. Nous avons induit ce que Pearson appelait une fausse corrélation (spurious correlation). En 1960, Chayes proposa que ces fausses corrélations soient induites non seulement sur des ratios de valeurs absolues, mais aussi sur des ratios d’une somme totale. Par exemple, dans une composition simple de deux types d’utilisation du territoire, si une proportion augmente, l’autre doit nécessairement diminuer.

n <- 30
tibble(A = runif(n, 0, 1)) %>%
  mutate(B = 1 - A) %>%
  ggplot(aes(x=A, y=B)) +
  geom_point() +
  coord_equal()
Complémentarité intrinsèque des données compositionnelles.

Figure 1.4: Complémentarité intrinsèque des données compositionnelles.

1.2 Les données compositionnelles

Les variables exprimées relativement à une somme totale sont dites compositionnelles. Les variables compositionnelles peuvent être des concentrations, des décomptes, des dimensions (aires ou volumes), des masses, des dépenses monétaires: en fait, tout ce qui s’exprime en proportions d’un tout. Elles possèdent les caractéristiques suivantes.

  1. Redondance d’information. Un système de deux proportions ne contient qu’une seule variable du fait que l’on puisse déduire l’une en soutrayant l’autre de la somme totale. Un vecteur compositionnel contient de l’information redondante. Pourtant, effectuer des statistiques sur l’une plutôt que sur l’autre donnera des résultats différents.
  2. Dépendance d’échelle. Les statistiques devraient être indépendantes de la somme totale utilisée. Pourtant, elles différeront sur l’on utilise par exemple, une proportion des mâles d’une part et des femelles d’autre part, ou la proportion de la somme des deux, de même que les résultats d’un test sanguin différera si l’on utilise une base sèche ou une base humide.
  3. Distribution théorique des données. Étant donnée que les proportions sont confinées entre 0 et 1 (ou 100%, ou une somme totale quelconque), la distribution normale (qui s’étend de -∞ à +∞) n’est souvent pas appropriée. On pourra utiliser la distribution de Dirichlet ou la distribution logitique-normale, mais d’autres approches sont souvent plus pratiques.

Pour illustrer l’effet de la distribution, voyons un diagramme ternaire incluant le sable, le limon et l’argile. En utilisant des écart-types univariés, nous obtenons l’ellipse en rouge, qui non seulement représente peu l’étalement des données, mais elle dépasse les bornes du triangle, admettant ainsi des proportions négatives. En bleu, la distribution logistique normale (issue des méthodes présentées plus loin dans cette section) convient davantage.

Représentation ternaire d'un écart-type multivarié biaisé (rouge) et non-biaisé (en bleu).

Figure 1.5: Représentation ternaire d’un écart-type multivarié biaisé (rouge) et non-biaisé (en bleu).

Les conséquences d’effectuer des statistiques linéaires sur des données compositionnelles brutes peuvent être majeures. En outre, Pawlowksy-Glahn et Egozcue (2006), s’appuyant en outre sur Rock (1988), note les problèmes suivants (exprimés en mes mots).

  1. les régressions, les regroupements et les analyses en composantes principales peuvent avoir peu ou pas de signification
  2. les propriétés des distributions peuvent être générées par l’opération de fermeture de la composition (s’assurer que le total des proportions donne 100%)
  3. les résultats d’analyses discriminantes linéaires sont propices à être illusoires
  4. tous les coefficients de corrélation seront affectés à des degrés inconnus
  5. les résultats des tests d’hypothèses seront intrinsèquement faussés

Pour contourner ces problèmes, il faut d’abord aborder les données compositionnelles pour ce qu’elles sont: des données intrinsèquement multivariées. Elles sont un nuage de point multidimensionnel, et non pas une collection de variables individuelles. Ceci qui n’empêche pas d’effectuer des analyses consciencieusement sous des angles particuliers.

Chargeons le module compositions (n’oubliez pas de l’installer au préalable) pour accéder à des données fictives de proportions de sable, limon et argile dans des sédiments. Notez que j’utiliserai la programmation en pipelines typique de tidyverse.

library("compositions")
data("ArcticLake")
ArcticLake <- ArcticLake %>% as_tibble()
head(ArcticLake)
## # A tibble: 6 x 4
##    sand  silt  clay depth
##   <dbl> <dbl> <dbl> <dbl>
## 1  77.5  19.5   3    10.4
## 2  71.9  24.9   3.2  11.7
## 3  50.7  36.1  13.2  12.8
## 4  52.2  40.9   6.6  13  
## 5  70    26.5   3.5  15.7
## 6  66.5  32.2   1.3  16.3

En R, on pourra aisément rapporter une composition en somme unitaire grâce à la fonction base::apply(). Pour information, la fonction base::apply() applique la fonction FUN = x/sum(x) pour toutes les lignes (MARGIN = 1, pour les colonnes, ce serait MARGIN = 2). Les résultats étant rendus sous forme de colonne, on doit finalement transposer (base:t()).

comp <- ArcticLake %>%
  dplyr::select(-depth) %>%
  apply(., MARGIN = 1, FUN = function(x) x/sum(x)) %>%
  t(.)

Plus simplement, pourra aussi utiliser la fonction compositions::acomp() (acomp pour Aitchison-composition) pour fermer la composition à une somme de 1.

comp <- ArcticLake %>%
  dplyr::select(-depth) %>%
  acomp(.)
comp[1:6, ]
##           sand      silt      clay
## [1,] 0.7750000 0.1950000 0.0300000
## [2,] 0.7190000 0.2490000 0.0320000
## [3,] 0.5070000 0.3610000 0.1320000
## [4,] 0.5235707 0.4102307 0.0661986
## [5,] 0.7000000 0.2650000 0.0350000
## [6,] 0.6650000 0.3220000 0.0130000

Cette stratégie a pour avantage d’attribuer à la variable comp la classe acomp, qui automatise les opérations dans l’espace compositionnel (que l’on nomme aussi le simplex). Notamment, la représentation ternaire est souvent utilisée pour présenter des compositions. Toutefois, il est difficile d’interpréter les compositions de plus de trois parties. La classe acomp automatise aussi la représentation teranaire.

plot(comp)
Représentation ternaire d'une composition granulométrique.

Figure 1.6: Représentation ternaire d’une composition granulométrique.

1.3 Les transformations compositionnelles

Afin de transposer cet espace clôt entre 0 et 1 en un espace ouvert de \(-\infty\) à \(+\infty\), on pourra diviser chaque proportion par une proportion de référence choisie parmi n’importe quelle proportion. Du coup, on retire une dimension redondante! En utilisant le log du ratio, l’inverse du ratio ne sera qu’un changement de signe, ce qui est pratique en statistiques linéaries. Cette solution, proposée par Aitchison (1986), est valide peu importe le nombre d’éléments formant la compositions. Pour une composition de \(A\), \(B\), \(C\), \(D\) et \(E\):

\[alr_A = log \left( \frac{A}{E} \right), alr_B = log \left( \frac{B}{E} \right), alr_C = log \left( \frac{C}{E} \right), alr_D = log \left( \frac{D}{E} \right)\]

Dans R, la colonne de référence est par défaut la dernière colonne de la matrice des compositions.

add_lr <- alr(comp)

Cette dernière stratégie se nomme les log-ratios additifs (\(alr\) pour additive log-ratio). Bien que valide pour effectuer des tests statistiques, cette stratégie a le désavantage de dépendre de la décision arbitraire de la composante à utiliser au numérateur. Deuxième restriction des alr: les axes de l’espace des alr n’étant pas orthogonaux, ils ne peuvent pas être utilisés pour effectuer des statistiques basées sur les distances euclidiennes étant donnée que la distance entre les points dépend de la variable choisie comme dénominateur commun.

add_lr_switch <- alr(comp[, c(2, 3, 1)])
dist(unclass(add_lr[1:6, ]), method = "euclidean")
##           1         2         3         4         5
## 2 0.2276858                                        
## 3 2.0933591 2.0527035                              
## 4 1.1846152 1.0686434 1.0912475                    
## 5 0.2979638 0.1195600 1.9389490 0.9491006          
## 6 1.5021428 1.4204061 3.3998698 2.3248142 1.5121727
dist(unclass(add_lr_switch[1:6, ]), method = "euclidean")
##           1         2         3         4         5
## 2 0.3486014                                        
## 3 2.1713475 1.9078113                              
## 4 1.6405267 1.3254326 0.7286134                    
## 5 0.4820630 0.1465561 1.7668173 1.1788856          
## 6 0.9461828 0.8883672 2.6176853 1.9278663 0.9708201

Les statistiques linéaires sont toutefois les mêmes peu importe le dénominateur commun: elles seront seulement exprimées selon les axes définis par la composition de référence (dénominateur commum).

Une stratégie proposée par Aitchison qui permet la mesure de distances euclidienne est d’effectuer un log-ratio entre chaque composante et la moyenne géométrique de toutes les composantes. Cette transformation se nomme le log-ratio centré (\(clr\), pour centered log-ratio)

\[clr_i = log \left( \frac{x_i}{g \left( x \right)} \right)\]

En R,

cen_lr <- clr(comp)

Avec des CLRs, les distances sont valides. Mais… nous restons avec le problème de la redondance d’information. En fait, la somme de chacunes des lignes d’une matrice de clr est de 0. Pas très pratique lorsque l’on effectue des statistiques incluant une inversion de la matrice de covariance (distance de Mahalanobis, géostatistiques, etc.)

cen_lr %>%
  cov() %>%
  solve()
 Error in solve.default(.) : le système est numériquement singulier : conditionnement de la réciproque = 4.44407e-17

Enfin, une autre méthode de transformation développée par Egoscue et al. (2003), les log-ratios isométriques (ou isometric log-ratios, ilr) projette les compositions comprenant D composantes dans un espace restreint de D-1 dimensions orthonormées. Ces dimensions doivent doivent être préalablement établie dans un dendrogramme de bifurcation, où chaque composante ou groupe de composante est successivement divisé en deux embranchement. La manière d’arranger ces balances importe peu, mais on aura avantage à créer des balances interprétables.

Le diagramme de balances, ou dendrogramme compositionnel, peut être encodé dans une partition binaire séquentielle (ou sequential bianry partition, sbp). Une sbp est une matrice de contrastes orthogonaux ou chaque ligne représente une partition entre deux variables ou groupes de variables. Une composante étiquettée +1 correspondra au groupe du numérateur, une composante étiquettée -1 au dénominateur et une composante étiquettée 0 sera exclue de la partition (Parent et al., 2013). J’ai reformulé la fonction CoDaDendrogram pour que l’on puisse ajouter des informations intéressantes sur les balants horizontaux. Cette fonction est disponible sur github.

source("https://raw.githubusercontent.com/essicolo/AgFun/master/codadend2.R")

sbp <- matrix(c(1, 1,-1,
                1,-1, 0),
              byrow = TRUE,
              ncol = 3)

CoDaDendrogram2(comp,
                V = gsi.buildilrBase(t(sbp)),
                ylim = c(0, 1),
                equal.height = TRUE,
                type = "boxplot",
                group = cut(ArcticLake$depth, 3))
Dendrogramme compositionnel: de la surface (rouge) vers la profondeur (bleu).

Figure 1.7: Dendrogramme compositionnel: de la surface (rouge) vers la profondeur (bleu).

Si la SBP est plus imposante (la précédente n’a que trois composantes, mais en pratique il y en a souvent davantage), il pourrait être plus aisé de monter dans un chiffrier, puis de l’importer dans R via un fichier csv.

Le calcul des ILRs est effectué comme suit.

\[ilr_j = \sqrt{\frac{n_j^+ n_j^-}{n_j^+ + n_j^-}} log \left( \frac{g \left( c_j^+ \right)}{g \left( c_j^+ \right)} \right)\]

où, à la ligne \(j\) de la SBP, \(n_j^+\) et \(n_j^-\) sont respectivement le nombre de composantes au numérateur et au dénominateur, \(g \left( c_j^+ \right)\) est la moyenne géométrique des composantes au numérateur et \(g \left( c_j^- \right)\) est la moyenne géométrique des composantes au dénominateur.

Les balances sont conventionnellement notées [A,B | C,D], ou les composantes A et B au dénominateur sont balancées avec les composantes C and D au numérateur. Une balance positive signifie que la moyenne géométrique des concentrations au numérateur est supérieur à celle au dénominateur, et inversement, alors qu’une balance nulle signifie que les moyennes géométriques sont égales (équilibre). Ainsi, en modélisation linéaire, un coefficient positif sur [A,B | C,D] signifie que l’augmentation de l’importance de C et D comparativement à A et B est associé à une augmentation de la variable réponse du modèle.

En R,

iso_lr <- ilr(comp, V = gsi.buildilrBase(t(sbp)))

Notez la forme gsi.buildilrBase(t(sbp)) est une opération pour obtenir la matrice d’orthonormalité à partir de la SBP.

Les ILRs sont des balances multivariées sur lesquelles on pourra effectuer des statistiques linéaries. Bien que l’interprétation des résultats comme collection d’interprétations sur des balances univariées pourra être affectée par la structure de la SBP, ni les statistiques linéaires multivariées, ni la distance entre les points ne seront affectés. En effet, chaque variante de la SBP est une rotation (d’un facteur de 60°) par rapport à l’origine:

source("lib/ilr-rotation-sbp.R")
Rotation des axes obtenus selon trois SBP différentes.

Figure 1.8: Rotation des axes obtenus selon trois SBP différentes.

En somme, la transformation par log-ratio consiste à exploser les données dans un espace libéré des contraintes compositionnelles.

C’est exactement ça, Bugs!

Pour les transformations inverses, c’est-à-dire pour passer de données transformées à des données exprimées dans l’échelle originale, vous pourrez utiliser les fonctions alrInv, clrInv et ilrInv. Dans tous les cas, si vous tenez à garder la trace de vos données dans leur format original, vous aurez avantage à ajouter à votre vecteur compositionnel la valeur de remplissage, constitué d’un amalgame des composantes non mesurées. Par exemple,

pourc <- c(N = 0.03, P = 0.001, K = 0.01)
acomp(pourc) # vous perdez la trace des proportions originales
##          N          P          K 
## 0.73170732 0.02439024 0.24390244 
## attr(,"class")
## [1] acomp
pourc <- c(N = 0.03, P = 0.001, K = 0.01)
Fv <- 1 - sum(pourc)
comp <- acomp(c(pourc, Fv = Fv))
comp
##     N     P     K    Fv 
## 0.030 0.001 0.010 0.959 
## attr(,"class")
## [1] acomp
iso_lr <- ilr(comp) # avec une sbp par défaut
ilrInv(iso_lr)
##      1    2     3    4    
## [1,] 0.03 0.001 0.01 0.959
## attr(,"class")
## [1] acomp

Les variables générée par les transformations compositionnelles doivent être interprétées avec précaution.

  • Un alr (log-ratio additif) est le log d’une proportion relative à une autre. Un vecteur d’alr est un point projeté dans un hyperespace compositionnel dont les axes sont à angle aigu (60 degrés): les distances entre les points ne sont pas euclidiennes.
  • Un clr (log-ratio centré) est le log d’une propostion relative à la moyenne (géométrique) des proportions. Les angles entre les axes d’un vecteur de clr sont orthogonaux. Néanmoins la somme de chaque vecteur de clr étant de 0, sa matrice de covariance ne peut pas être inversée (singularité), ce qui rend notamment impossible le calcul de la distance de Mahalanobis et de certains autres calculs géostatistiques.
  • Un ilr (log-ratio isométrique) est le log normé du ration entre deux proportions ou moyenne géométrique d’un ensemble de proportions. Les angles entre les axes d’un vecteur de ilr sont orthogonaux et la matrice de covariance peut être inversée, mais la définition des axes comme balance entre des proportions ou groupes de proportions est souvent difficile à interpréter.

Étant donnée que les données compositionnelles sont intrinsèquement multivariées, chaque variable qui compose un vecteur compositionnel, transformé ou non, constitue une perspective particulière d’un nuage de points. Or une collection des perspectives individuelles reste une appréciation grossière d’un espace multivarié.

La transformation à privilégié dépendra de l’analyse que vous entreprendrez. Toutefois, Filzmoser et al. (2018) recommande les transformations ilr.

[The] Euclidean vector space structure [of isometric log-ratios] enables to construct coordinates with respect to a basis, eventually coefficients of a generating system. Here, isometric logratio coordinates, real coordinates with respect to an orthonormal basis in the Aitchison geometry, are preferable. - Filzmoser et al. (2018)

1.4 Les zéros

Les proportions négatives ou plus imposantes que le tout sont exclues. De même, une proportion de 100%, occupant tout le simplex, implique que les proportions complémentaires soient nulles. Or, le log de 0 est \(-\infty\), ce qui dispose de la possibilité d’effectuer des transformations logarithmiques. En fait, les données compositionnelles sont exprimées entre 0 et la somme totale, ce qui exclut la possibilité d’obtenir des proportions de 0% ou de 100%. Or, il arrive souvent qu’un tableau de données inclut des 0. Dans ce cas, on doit mener une réflexion sur la nature des valeurs nulles.

1.4.1 Les zéros arrondis

Arrondissez 0.01 à une seule décimale et vous obtiendrez 0.0. Ou bien, mesurez une concentration de 0.01 avec un appareil dont le seuil de détection est 0.1, vous obtiendrez encore 0.0. Étant des valeurs qui ne peuvent être observées pour des raisons méthodologiques, les zéros arrondis sont équivalents à des valeurs manquantes non aléatoires (Martin-Fernandez et al., 2011). Il existe plusieurs manières de remplacer ces zéros par une valeur attendue. La méthode la plus simple est de remplacer les zéros par une proportion du seuil de détection, par exemple 65% (Martin-Fernandez et al., 2003).

Chargeons d’abord des données fictives.

library("zCompositions")
set.seed(542353)
data(LPdata)
rows <- sample(1:nrow(LPdata), 3)
LPdata[rows, ]
##      Cr  B   P  V  Cu   Ti   Ni  Y Sr La  Ce Ba Li     K Rb
## 38 27.5 17 148 29 2.7 4335  0.0 17 24  8 163 74 30  9703 48
## 3  30.4 23 433 42 3.8 3305 16.6 22 59 14 240 75 80 12209 53
## 79 25.6 14 135 33 0.0 3925 14.2 20 18  9  97 36 77  3319 41

Si le seuil de détection est de 1 pour toutes les colonnes,

LPdata_z <- LPdata
LPdata_z[LPdata == 0] <- 1*0.65
LPdata_z[rows, ]
##      Cr  B   P  V   Cu   Ti    Ni  Y Sr La  Ce Ba Li     K Rb
## 38 27.5 17 148 29 2.70 4335  0.65 17 24  8 163 74 30  9703 48
## 3  30.4 23 433 42 3.80 3305 16.60 22 59 14 240 75 80 12209 53
## 79 25.6 14 135 33 0.65 3925 14.20 20 18  9  97 36 77  3319 41

Par exemple, le remplacement multiplicatif simple incluse dans le module zCompositions, dont la mathématique varie légèrement.

LPdata_z <- multRepl(LPdata, dl = rep(1, ncol(LPdata)), delta = 0.65, label = 0)
LPdata_z[rows, ]
##      Cr  B   P  V        Cu   Ti         Ni  Y Sr La  Ce Ba Li     K Rb
## 38 27.5 17 148 29 2.7000000 4335  0.6500289 17 24  8 163 74 30  9703 48
## 3  30.4 23 433 42 3.8000000 3305 16.6000000 22 59 14 240 75 80 12209 53
## 79 25.6 14 135 33 0.6500544 3925 14.2000000 20 18  9  97 36 77  3319 41

1.4.2 Les décomptes de zéro

Vous décomptez des espèces sur une parcelle ou bien séquencez des OTU (operational taxonomy units, qui peuvent être interprétés - avec une précaution bioinformaticienne - comme des espèces): certaines espèces ne seront pas observées, bien que leur occurence puisse être attendue. Si vous vous intéressez davantage aux proportions d’espèces qu’à leur décompte absolu, vous devriez transformer vos données avec les technique présentées ci-dessus. De même que précédemment, la non observation d’un décompte peut provienir de contraintes méthodologiques: parcelle trop petite, temps d’observation trop court, ou quoi qui pourrait incluencer la non observation d’une fréquence rare. Par exemple, les décomptes de cochons occupés à différentes activités.

set.seed(30147)
data(Pigs)
rows <- sample(1:nrow(Pigs), 3)
Pigs[rows, ]
##    BED HALF.BED PASSAGE HALF.PASS FEEDER HALF.FEED
## 19  83        1       9         1      3         0
## 15  68        1      16         1     10         1
## 18  69        0      26         0      2         0

En ces circonstances, il est possible de remplacer les zéros par des méthodes bayésiennes. Les données sont du coup transformées en proportions.

Pigs_z <- cmultRepl(Pigs, method = "CZM") # CZM: count zero multiplicative
Pigs_z[rows, ]
##          BED    HALF.BED    PASSAGE   HALF.PASS     FEEDER   HALF.FEED
## 19 0.8528032 0.010274737 0.09247263 0.010274737 0.03082421 0.003350515
## 15 0.7010309 0.010309278 0.16494845 0.010309278 0.10309278 0.010309278
## 18 0.7041901 0.003350515 0.26534701 0.003350515 0.02041131 0.003350515

Dans le cas des OTU, il arrive souvent d’obtenir au-dessus de 90% de zéros. Une option est de regrouper (amalgamer) autant que possible les OTU en taxons d’ordre supérieurs.

1.4.3 Les zéros structurels

Il peut arriver que les zéros soient générés par les mécanismes étudiés, non pas des erreurs de mesure. Par exemple, certaines espèces peuvent ne pas occuper certaines aires cadrillées parce qu’elles sont structurellement exclues de leur niche écologique (une épinette dans un lac ou la proportion de viande au menu d’une famille véganne). Pour ce type de zéro, il n’y a pas de réponse facile. De même que précédemment, il pourrait être possible de procéder à l’amalgamation. Une méthode conviviale est de séparer l’analyse en différents groupes (par exemple, on décompte généralement des espèces différentes sur terre et dans un lac).

Si vos données font partie d’un tout, je vous recommande chaudement d’utiliser des méthodes compositionnelles autant pour l’analyse que la modélisation. Parmi les trois transformations présentées, Pour en savoir davantage, le livre Compositional data analysis with R, de van den Boogart et Tolosana-Delgado (2013) pourra vous guider, quoi que je recommande davantage le plus récent et à mon avis mieux vulgarisé [Applied Compositional Data Analysis With Worked Examples in R], de Filzmoser et al. (2018).