5.2 Distance entre individus

Pour partir d’un exemple concret, imaginez que vous êtes en train d’analyser des données concernant les échantillons de plancton que vous avez prélevé sur votre lieu de recherche. Ce plancton a été numérisé (photo de chaque organisme) et les images ont été traitées avec un logiciel qui mesure automatiquement une vingtaine de variables telles que la surface de l’objet sur l’image, son périmètre, sa longueur, … Vous vous trouvez donc face à un jeu de données qui a une taille non négligeable : 20 colonnes par 1262 lignes, soit le nombre d’individus mesurés dans vos échantillons.

zoo <- read("zooplankton", package = "data.io")
zoo

# # A tibble: 1,262 x 20
#      ecd  area perimeter feret major minor  mean  mode   min   max std_dev range
#    <dbl> <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>
#  1 0.770 0.465      4.45 1.32  1.16  0.509 0.363 0.036 0.004 0.908   0.231 0.904
#  2 0.700 0.385      2.32 0.728 0.713 0.688 0.361 0.492 0.024 0.676   0.183 0.652
#  3 0.815 0.521      4.15 1.33  1.11  0.598 0.308 0.032 0.008 0.696   0.204 0.688
#  4 0.785 0.484      4.44 1.78  1.56  0.394 0.332 0.036 0.004 0.728   0.218 0.724
#  5 0.361 0.103      1.71 0.739 0.694 0.188 0.153 0.016 0.008 0.452   0.110 0.444
#  6 0.832 0.544      5.27 1.66  1.36  0.511 0.371 0.02  0.004 0.844   0.268 0.84 
#  7 1.23  1.20      15.7  3.92  1.37  1.11  0.217 0.012 0.004 0.784   0.214 0.78 
#  8 0.620 0.302      3.98 1.19  1.04  0.370 0.316 0.012 0.004 0.756   0.246 0.752
#  9 1.19  1.12      15.3  3.85  1.34  1.06  0.176 0.012 0.004 0.728   0.172 0.724
# 10 1.04  0.856      7.60 1.89  1.66  0.656 0.404 0.044 0.004 0.88    0.264 0.876
# # … with 1,252 more rows, and 8 more variables: size <dbl>, aspect <dbl>,
# #   elongation <dbl>, compactness <dbl>, transparency <dbl>, circularity <dbl>,
# #   density <dbl>, class <fct>

Vous voulez regrouper votre plancton en fonction de la ressemblance entre les organismes, c’est-à-dire, en fonction des écarts entre les mesures effectuées pour les 19 variables quantitatives, à l’exclusion de la vingtième colonne class qui est une variable factor). En raison de la taille du tableau, il est évident que cela ne pourra pas se faire de manière manuelle. Nous pouvons raisonnablement considérer que plus les mesures sont similaires entre deux individus, plus ils ont des chances d’être semblables, c’est-à-dire, d’appartenir au même groupe taxonomique. Mais comment faire pour synthétiser l’information de similarité ou différence contenue dans 19 paires de valeurs (une paire par variable)  ? Nous avons besoin d’une mesure de distance qui quantifie la similarité (ou à l’inverse la dissimilarité) en un seul nombre. Celle qui vient naturellement à l’esprit est la distance euclidienne. Prenons un cas simplifié. Quelle est la distance qui sépare deux individus A et B par rapport à trois variables x, y, z ? Ici, nous pouvons représenter l’information graphiquement dans un espace à trois dimensions. La distance qui nous intéresse est la distance linéaire entre les deux points dans l’espace. Autrement dit, c’est la longueur du segment de droite qui relie les deux points dans l’espace. Cette distance, nous pouvons la calculer à l’aide de la formule suivante (voir par exemple ici pour une résolution dans le plan) :

\[\mathrm{D_{Euclidean}}_{A, B} = \sqrt{(x_A - x_B)^2 + (y_A - y_B)^2 + (z_A - z_B)^2}\]

Notez que cette formule se généralise à n dimensions et s’écrit alors, pour n’importe quelle paire d’individus indicés j et k dans notre tableau et pour les différentes mesures de 1 à i notées y_i :

\[\mathrm{D_{Euclidean}}_{j, k} = \sqrt{\sum_{i=1}^{n}(y_{ij}-y_{ik})^2}\]

C’est la racine carré de la somme dans les n dimensions des écarts entre les valeurs au carré pour toutes les variables y_i. Plus sa valeur est grande, plus les individus sont éloignés (différents). Pour cette raison, nous appellerons cette distance, une mesure de dissimilarité.

5.2.1 Matrice de distances

Nous avons maintenant la possibilité de quantifier la similitude entre nos organismes planctoniques… mais nous en avons un grand nombre. Cela va être impossible à gérer autant de mesures qu’il y a de paires possibles parmi 1262 individus¹². La matrice de distance est une matrice ici 1262 par 1262 qui rassemble toutes les valeurs possibles. Notez que sur la diagonale, nous comparons chaque individu avec lui-même. La distance euclidienne vaut donc systématiquement zéro sur la diagonale.

\[\mathrm{D_{Euclidean}}_{j, j} = 0\]

De plus, de part et d’autre de cette diagonale, nous trouvons les paires complémentaires (j versus k d’un côté et k versus j de l’autre). Or qu’elle soit mesurée dans un sens ou dans l’autre, la distance du segment de droite qui relie deux points dans l’espace est toujours la même.

\[\mathrm{D_{Euclidean}}_{j, k} = \mathrm{D_{Euclidean}}_{k, j}\]

Par conséquent, seulement une portion (soit le triangle inférieur, soit le triangle supérieur hors diagonale) est informative. La diagonale ne porte aucune information utile, et l’autre triangle est redondant. Nous avons donc pour habitude de ne calculer et représenter que le triangle inférieur de cette matrice.

À vous de jouer !

Avant de poursuivre, nous allons utiliser des fonctions personnalisées qui nous faciliteront le travail. Le code qui suit (dépliez la zone de code si vous êtes curieux) va nous créer ces fonctions. Pour un document R Markdown, nous mettrons ce code dans un script R séparé, placé dans R/functions.R dans notre projet, et nous ferons source("../R/functions.R") dans un chunk de setup. A terme, ces fonctions seront intégrées directement dans la surcouche SciViews::R par dessus R.

# CAH for SciViews, version 1.2.0
# Copyright (c) 2021, Philippe Grosjean (phgrosjean@sciviews.org)

SciViews::R

# dist is really a dissimilarity matrix => we use dissimilarity() as in the
# {cluster} package, i.e., class is c("dissimilarity", "dist")
# TODO: also make a similarity object and convert between the two
# fun can be stats::dist, vegan::vegdist, vegan::designdist, cluster::daisy
# factoextra::get_dist and probably other dist-compatible functions
# Depending on method =, use either vegan::vegdist or stats::dist as default fun
dissimilarity <- function(data, formula = ~ ., subset = NULL,
  method = "euclidean", scale = FALSE, rownames.col = "rowname",
  transpose = FALSE, fun = NULL, ...) {
  # TODO: get more meaningful warnings and errors by replacing fun by actual
  # name of the function
  if (is.null(fun)) {# Default function depends on the chosen method
    if (method %in% c("maximum", "binary", "minkowski")) {
      fun <- stats::dist
    } else {
      fun <- vegan::vegdist # Note: euclidean/manhattan/canberra in both, but
      # we prioritize vegdist, and canberra is not calculated the same in dist!
    }
  }
  # We accept only formulas with right-hand side => length must be two
  if (length(formula) == 3)
    stop("The formula cannot have a left-hand term")

  # With matrices, we don't use rownames.col: rownames are already correctly set
  if (!is.matrix(data)) {# row names may be in a column (usual for tibbles)
    data <- as.data.frame(data)
    if (rownames.col %in% names(data)) {
      rownames(data) <- data[[rownames.col]]
      data[[rownames.col]] <- NULL
    } else {# rownames.col is NOT used
      rownames.col <- NULL
    }
    if (as.character(formula[2] != ".")) {
      # Subset the columns
      data <- model.frame(formula, data = data, subset = subset)
    } else if (!is.null(subset)) {
      data <- data[subset, ]
    }
  } else {# A matrix
    rownames.col <- NULL
    if (as.character(formula[2] != ".")) {
      # Subset the columns (and possibly the rows)
      if (is.null(subset)) {
        data <- data[, all.vars(formula)]
      } else {
        data <- data[subset, all.vars(formula)]
      }
    }
  }

  if (isTRUE(transpose))
    data <- t(data)

  # Arguments method =/metric = and stand = not always there
  if (!is.null(as.list(args(fun))$metric)) {# metric = instead of method =
    dst <- fun(data, metric = method, stand = scale, ...)
  } else if (isTRUE(scale)) {
    if (is.null(as.list(args(fun))$stand)) {# fun has no stand = argument
      data <- scale(data)
      dst <- fun(data, method = method, ...)
    } else {# We don't standardise ourself because there may be also qualitative
      # or binary data (like for cluster::daisy, for instance)
      dst <- fun(data, method = method, stand = scale, ...)
    }
  } else {# Just method = and scale = FALSE
    dst <- fun(data, method = method, ...)
  }
  attr(dst, "call") <- match.call()
  # Depending if it is a dist or dissimilarity object, the method is stored in
  # method or in Metric, but we use metric in our own version to avoid a clash
  # with the method item in cluster()/hclust() further on (hclust change it
  # into dist.method, but it is better to have the right name right now)
  attr(dst, "metric") <- method
  # dist or dissimilarity object use Labels, but we use labels everywhere else
  # including in cluster()/hclust()
  # So, we make sure labels is present (in hclust, it is labels anyway!)
  attr(dst, "labels") <- rownames(data)
  # Default values for Diag and Upper set to FALSE
  if (is.null(attr(dst, "Diag"))) attr(dst, "Diag") <- FALSE
  if (is.null(attr(dst, "Upper"))) attr(dst, "Upper") <- FALSE
  # Keep info about how raw data were transformed
  attr(dst, "rownames.col") <- rownames.col
  attr(dst, "transpose") <- transpose
  attr(dst, "scale") <- scale
  class(dst) <- unique(c("dissimilarity", class(dst)))
  dst
}

as.dissimilarity <- function(x, ...)
  UseMethod("as.dissimilarity")
as_dissimilarity <- as.dissimilarity # Synonym

as.dissimilarity.matrix <- function(x, ...) {
  dst <- as.dist(x, ...)
  attr(dst, "call") <- match.call()
  attr(dst, "metric") <- attr(dst, "method") # Make sur metric is used
  class(dst) <- unique(c("dissimilarity", class(dst)))
  dst
}

# We want to print only the first few rows and columns
print.dissimilarity <- function(x, digits.d = 3L, rownames.lab = "labels",
...) {
  mat <- as.matrix(x)
  mat <- format(round(mat, digits.d))
  diag(mat) <- ""
  mat[upper.tri(mat)] <- ""
  class(mat) <- c("dst", "matrix")
  tbl <- tibble::as_tibble(mat)
  #tbl <- tibble::add_column(tbl, {{rownames.lab}} = rownames(mat), .before = 1)
  # I prefer this
  tbl <- dplyr::bind_cols(
    as_tibble_col(rownames(mat), column_name = rownames.lab), tbl)
  tbl <- tbl[, -ncol(tbl)]
  more_info <- ""
  if (isTRUE(attr(x, "scale"))) {
    if (isTRUE(attr(x, "transpose"))) {
      more_info <- " (transposed then scaled data)"
    } else {# Only scaled
      more_info <- " (scaled data)"
    }
  } else {
    if (isTRUE(attr(x, "transpose")))
      more_info <- " (transposed data)"
  }
  cat("Dissimilarity matrix with metric: ", attr(x, "metric"),
    more_info, "\n", sep = "")
  print(tbl)
  invisible(x)
}

labels.dissimilarity <- function(object, ...) {
  labs <- object$labels
  if (is.null(labs)) object$Labels
}

nobs.dissimilarity <- function(object, ...)
  attr(object, "Size")

# TODO: `[` by first transforming into a matrix with as.matrix()

autoplot.dissimilarity <- function(object, order = TRUE, show.labels = TRUE,
lab.size = NULL, gradient = list(low = "red", mid = "white", high = "blue"),
...) {
  factoextra::fviz_dist(object, order = order, show_labels = show.labels,
    lab_size = lab.size, gradient = gradient)
}

chart.dissimilarity <- function(data, ...,
type = NULL, env = parent.frame())
  autoplot(data, type = type, ...)

# cluster object (inheriting from hclust)
cluster <- function(x, ...)
  UseMethod("cluster")

cluster.default <- function(x, ...)
  stop("No method for object of class ", class(x)[1])

# Cluster uses hclust() by default, ... but it looks first for a faster
# implementation in either {fastcluster} or {flashClust} before falling back
# to the {stats} version.
# The functions cluster::agnes() and cluster::diana() should be compatible too,
# as well as any function that returns an object convertible into hclust
# by as.hclust() (but not tested yet)
# Also, a version where the raw data are provided and the disimilarity matrix
# is internally calculated should be also provided (see cluster::agnes)
# See also {ape} for phylogenetic trees methods
cluster.dist <- function(x, method = "complete", fun = NULL, ...) {
  if (is.null(fun)) {
    # We try fastcluster, then flashClust, then stats
    fun <- try(fastcluster::hclust, silent = TRUE)
    if (inherits(fun, "try-error"))
      fun <- try(flashClust::hclust, silent = TRUE)
    if (inherits(fun, "try-error"))
      fun <- try(stats::hclust, silent = TRUE)
  }
  clst <- fun(x, method = method, ...)
  clst <- as.hclust(clst)
  clst$call <- match.call()
  # hclust has to give a different name to the distance metric: dist.method
  # but we use metric. Again, keep both for maximum compatibility
  clst$metric <- clst$dist.method
  # If the original data were scaled or transposed, get the info also
  clst$rownames.col <- attr(x, "rownames.col")
  clst$scale <- attr(x, "scale")
  clst$transpose <- attr(x, "transpose")
  class(clst) <- unique(c("cluster", class(clst)))
  clst
}

# A couple of useful methods for our cluster object
# str() method is gathered from a dendrogram object
str.cluster <- function(object, max.level = NA, digits.d = 3L, ...)
  str(as.dendrogram(object), max.level = max.level, digits.d = digits.d, ...)

labels.cluster <- function(object, ...)
  object$labels

nobs.cluster <- function(object, ...)
  length(object$order)

# Other methods by first transforming into dendrogram: rev, reorder, order, [[

# cutree() is an explicit name, but it does not follow the rule of using
# known methods... and here, it really something that predict() is made for,
# except it cannot handle newdata =, but that argument is not in its definition
predict.cluster <- function(object, k = NULL, h = NULL, ...)
  cutree(object, k = k, h = h)

# There is no broom::glance() or broom::tidy() yet (what to put in it?),
# but broom:augment() should be nice = add the clusters as .fitted in the tibble
library(broom)
augment.cluster <- function(x, data, k = NULL, h = NULL, ...) {
  # Should we transpose the data (note: this is against augment() rules, but...)
  if (isTRUE(x$transpose)) {
    # We first have to make sure rownames are correct before the transposition
    if (!is.matrix(data) && !is.null(data[[x$rownames.col]])) {
      rownames(data) <- data[[x$rownames.col]]
      data[[x$rownames.col]] <- NULL
    }
    data <- t(data)
    msg <- "transposed data"
  } else {
    msg <- "data"
  }
  data <- as_tibble(data)

  # Get clusters
  clst <- predict(x, k = k, h = h, ...)
  if (nrow(data) != length(clst)) {
    stop("Different number of items in ", msg, " (",nrow(data) ,
      ") and in the clusters (", length(clst), ")")
  }
  tibble::add_column(data, .fitted = clst)
}

# Instead of the default plot.hclust(), we prefer the plot.dendrogram() version
# that allows for more and better variations of the dendrogram (horizontal or
# circular), see http://www.sthda.com/english/wiki
# /beautiful-dendrogram-visualizations-in-r-5-must-known-methods
# -unsupervised-machine-learning
plot.cluster <- function(x, y, labels = TRUE, hang = -1, check = TRUE,
type = "vertical", lab = "Height", ...) {
  type <- match.arg(type[1], c("vertical", "horizontal", "circular"))
  # type == "circular" is special because we need to transform as ape::phylo
  if (type == "circular") {
    if (!missing(hang))
      warning("'hang' is not used with a circular dendrogram")
    phylo <- ape::as.phylo(x)
    plot(phylo, type = "fan", font = 1, show.tip.label = labels, ...)
  } else {# Use plot.dendrogram() instead
    # We first convert into dendrogram objet, then we plot it
    # (better that plot.hclust())
    if (isTRUE(labels)) leaflab <- "perpendicular" else leaflab <- "none"
    dendro <- as.dendrogram(x, hang = hang, check = check)
    if (type == "horizontal") {
      plot(dendro, horiz = TRUE, leaflab = leaflab, xlab = lab, ...)
    } else {
      plot(dendro, horiz = FALSE, leaflab = leaflab, ylab = lab, ...)
    }
  }
}

# This is to draw circles in a plot (where to cut in a circular dendrogram)
# TODO: should be nice to do similar function for other symbols too in SciViews
circle <- function(x = 0, y = 0, d = 1, col = 0, lwd = 1, lty = 1, ...)
  symbols(x = x, y = y, circles = d / 2, fg = col, lwd = lwd, lty = lty,
    inches = FALSE, add = TRUE, ...)

# TODO: make sure the dendrogram is correct with different ggplot themes
autoplot.cluster <- function(object, labels = TRUE, type = "vertical",
circ.text.size = 3, theme = theme_sciviews(), xlab = "", ylab = "Height", ...) {
  if (is.null(type))
    type <- "vertical"
  type <- match.arg(type[1], c("vertical", "horizontal", "circular"))

  # Create the dendrogram
  ddata <- ggdendro::dendro_data(object, type = "rectangle")
  dendro <- ggplot(ggdendro::segment(ddata)) +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
    theme + xlab(xlab) + ylab(ylab)

  if (type == "circular") {
    if (isTRUE(labels)) {
      # Get labels (need one more to avoid last = first!)
      label_df <- tibble::tibble(labels = c(labels(object)[object$order], ""))
      xmax <- nobs(object) + 1
      label_df$id <- 1:xmax
      angle <-  360 * (label_df$id - 0.5) / xmax
      # Left or right?
      label_df$hjust <- ifelse(angle < 270 & angle > 90, 1, 0)
      # Angle for more readable text
      label_df$angle <- ifelse(angle < 270 & angle > 90, angle + 180, angle)
    }

    # Make the dendrogram circular
    dendro <- dendro +
      scale_x_reverse() +
      scale_y_reverse() +
      coord_polar(start = pi/2)
    if (isTRUE(labels))
      dendro <- dendro +
        geom_text(data = label_df,
          aes(x = id, y = -0.02, label = labels, hjust = hjust),
          size = circ.text.size, angle = label_df$angle, inherit.aes = FALSE)
    dendro <- dendro +
      theme(panel.border = element_blank(),
        axis.text = element_blank(),
        axis.line = element_blank(),
        axis.ticks.y = element_blank()) +
      ylab("")

  } else if (type == "vertical") {# Vertical dendrogram
    dendro <- dendro +
      scale_x_continuous(breaks = seq_along(ddata$labels$label),
        labels = ddata$labels$label) +
      scale_y_continuous(expand = expansion(mult = c(0, 0.02))) +
      theme(panel.border = element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.line.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(angle = 90, hjust = 0.5))
    if (!isTRUE(labels))
      dendro <- dendro +
        theme(axis.text.x = element_blank())

  } else {# Horizontal dendrogram
    dendro <- dendro +
      scale_x_continuous(breaks = seq_along(ddata$labels$label),
        labels = ddata$labels$label, position = "top") +
      scale_y_reverse(expand = expansion(mult = c(0.05, 0))) +
      coord_flip() +
      theme(panel.border = element_blank(),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank())
    if (!isTRUE(labels))
      dendro <- dendro +
        theme(axis.text.y = element_blank())
  }
  dendro
}

chart.cluster <- function(data, ...,
  type = NULL, env = parent.frame())
  autoplot(data, type = type, ...)

# To indicate where to cut in the dendrogram, one could use `geom_hline()`,
# but when the dendrogram is horizontal or circular, this is suprizing. So,
# I define geom_dendroline(h = ....)
geom_dendroline <- function(h, ...)
  geom_hline(yintercept = h, ...)

Voici les fonctions à disposition (nous verrons leur usage progressivement dans ce module) :

• dissimilarity(data, method = "euclidean", scale = FALSE, transpose = FALSE) : matrice de dissimilarité, voir ?vegan::vegdist
  • print(), chart(), labels() et nobs() sont disponibles
• cluster(x, method = "complete") : CAH à partir d’un objet “dissimilarity”, voir ?stats:hclust
  • print(), str(), plot(), chart(), labels(), nobs(), predict() et augment()
• chart(cluster) : visualise un cluster (dendrogramme). Utiliser + geom_dendroline(h = XX, color = "red") pour y visualiser la coupure
• predict(cluster, h = 5) ou predict(cluster, k = 5) extrait les groupes
• augment(data = df, cluster, h =|k = ) ajoute les groupes dans le tableau df

Voici comment nous calculons une matrice de distances zoo6_dist, ici sur un petit sous-ensemble de six lignes de notre jeu de données (nous devons aussi éliminer la colonne class qui ne contient pas de données numériques et qui ne nous intéresse pas pour le moment) :

zoo %>.%
  select(., -class) %>.% # Élimination de la colonne class
  head(., n = 6) -> zoo6 # Récupération des 6 premiers individus
zoo6_dist <- dissimilarity(zoo6, method = "euclidean")
zoo6_dist

# Dissimilarities :
#            1          2          3          4          5
# 2  8.2185826                                            
# 3  2.5649705  5.7320911                                 
# 4  0.8582142  7.8813537  2.2220079                      
# 5  4.8478629  4.2950067  3.0119468  4.6148173           
# 6  2.4269520 10.6197317  4.9228255  2.8477882  7.1766853
# 
# Metric :  
# Number of objects : 6

Nous voyons bien ici que R n’imprime que le triangle inférieur de notre matrice 6 par 6. Notez aussi que les objets dist de tailles plus réalistes que vous utiliserez dans vos analyses ne sont prévue pour être imprimées et visualisées telles quelles. Il s’agit seulement de la première étape vers une représentation utile qui sera réalisée à la page suivante, à l’aide de la classification hiérarchisée.

Nous verrons plus loin comment nous pouvons utiliser l’information que cette matrice de distances contient pour regrouper les individus de manière pertinente. Mais avant cela, nous avons besoin d’un peu de théorie pour bien comprendre quelle métrique choisir pour calculer nos distances et pourquoi. On parle aussi d’indices de similarité ou dissimilarité.

Attention : nous n’avons pas considéré ici les unités respectives de nos variables. Une surface (mm²) ou une longeur (mm) ne sont pas mesurées dans les mêmes unités. Nous risquons alors de donner plus de poids dans nos calculs aux variables qui présentent des valeurs élevées. Nous aurions le même effet si nous décidions par exemple d’exprimer une mesure longitudinale en µm au lieu de l’exprimer en mm. Dans ce cas, il vaut mieux standardiser d’abord le tableau (moyenne de zéro et écart type de un) selon les colonnes avant d’effectuer le calcul. Une illustration de cette approche sera discutée plus loin.

À vous de jouer !

5.2.2 Indices de (dis)similarité

Un indice de similarité (similarity index en anglais) est une descripteur statistique (nombre unique) de la similitude de deux échantillons ou individus représentés par plusieurs variables dans un échantillon multivarié. Un indice de similarité prend une valeur comprise entre 0 (différence totale) et 1 ou 100% (similitude totale). Un indice de dissimilarité} est le complément d’un indice de similarité (dis = 1 – sim) ; sa valeur est comprise entre 100% (différence totale) et 0 (similitude totale).

Attention : dans certains cas, un indice de dissimilarité peut varier de 0 à +\(\infty\)**. Il n’existe alors pas d’indice de similarité complémentaire. C’est le cas précisément de la distance euclidienne que nous avons exploré jusqu’ici.

Tous les indices de similarité / dissimilarité peuvent servir à construire des matrices de distances.

5.2.2.1 Indice de Bray-Curtis

L’indice de dissimilarité de Bray-Curtis, aussi appelé coefficient de Czecanowski est calculé comme suit :

\[\mathrm{D_{Bray-Curtis}}_{j,k}=\frac{\sum_{i=1}^{n}\left|y_{ij}-y_{ik}\right|}{\sum_{i=1}^{n}(y_{ij}+y_{ik})}\]

Dans SciViews R nous utiliserons dissimilarity(DF, method = "bray"). cet indice s’utilise pour mesurer, entre autres, la similitude entre échantillon sur base du dénombrement d’espèces. Si le nombre d’espèces est très variable (espèces dominantes versus espèces rares), nous devons transformer les données pour éviter de donner trop de poids aux espèces les plus abondantes (ex: \(log(x+1)\), double racine carrée, …).

Une caractéristique essentielle de cet indice (contrairement à la distance euclidienne) est que toute double absence n’est pas prise en compte dans le calcul. C’est souvent pertinent dans le cadre de son utilisation comme le dénombrement d’espèces. En effet, quelle information utile retire-t-on de doubles zéros dans un tableau répertoriant la faune belge pour le crocodile du Nil et le tigre de Sibérie par exemple ? Aucune ! Ils sont tous deux systématiquement absents des dénombrements, mais cette double absence n’apporte aucune information utile pour caractériser la faune belge par ailleurs.

L’indices de similarité de Bray-Curtis (sim) est complémentaire à l’indices de dssimilarité correspondant (dis tel que calculé ci-dessus) :

\[sim = 1 – dis\]

5.2.2.2 Indice de Canberra

L’indice de dissimilarité de Canberra est apparenté à l’indice de Bray-Curtis mais il pondère les individus en fonction du nombre d’occurrences afin de donner le même poids à chacune. Il se calcule comme suit :

\[\mathrm{D_{Canberra}}_{j,k}=\frac{1}{nz}\sum_{i'=1}^{nz}\frac{\left|y_{i'j}-y_{i'k}\right|}{\left|y_{i'j}\right|+\left|y_{i'k}\right|}\]

où \(nz\) est le nombre de valeurs non nulles simultanément dans le tableau de départ. Toutes les cas contribuent ici de manière égale. C’est un point positif, mais il faut faire attention à ce que cet indice a souvent tendance à donner, au contraire, trop d’importance aux dénombrements très rares observés une seule fois ou un petit nombre de fois ! Dans SciViews R, nous utiliserons dissimilarity(DF, method = "canberra").

Toute double absence n’est pas prise en compte ici également. Seuls les indices ne dépendant pas des doubles zéros sont utilisables pour des dénombrements d’espèces ou des présence-absence. Ainsi pour ce type de données, notre choix se portera sur :

• Bray-Curtis si l’on souhaite que le résultat soit dominé par les espèces les plus abondantes.

• Canberra si notre souhait est de donner la même importance à toutes les espèces, mais avec un risque d’importance exagérée des espèces rares par rapport au nombre relatif d’observations pour l’ensemble du jeu de données.

• Bray-Curtis sur données transformées (\(log(x+1)\) ou double racine carrée) pour un compromis entre les deux avec prise en compte de toutes les espèces, mais dépondération partielle des espèces les plus abondantes. C’est souvent un bon compromis.

Attention : Si les volumes échantillonnés entre stations ne sont pas comparables, il faut standardiser (moyenne nulle et écart type un) les données selon les échantillons avant de faire les calculs de distances. L’argument scale = TRUE pourra être ajouté à l’appel de dissimilarity() pour ce faire.

De même que pour Bray-Curtis, l’indice de similarité sim se calcule à partir de l’indice de dissimilarité dis tel que ci-dessus comme \(sim = 1 - dis\).

5.2.2.3 Distance Euclidienne

Nous savons déjà que c’est la distance géométrique entre les points dans un espace à n dimensions :

\[\mathrm{D_{Euclidean}}_{j,k}=\sqrt{\sum_{i=1}^{n}(y_{ij}-y_{ik})^2}\]

Dans SciViews R, cette distance peut être calculée avec dissimilarity(DF, method = "euclidean"). Attention : en anglais, c’est euclidean avec un “e”, pas euclidian ! Cet indice de dissimilarité est utile pour des mesures quantitatives, pour des données environnementales, etc. Il faut que les mesures soient toutes effectuées dans les mêmes unités. Si ce n’est pas le cas, nous devons les standardiser avant le calcul, avec scale = TRUE. Il n’existe pas d’indice de similarité complémentaire.

5.2.2.4 Distance de Manhattan

La distance de Manhattan, encore appelée “city-block distance” en anglais, est un indice de dissimilarité qui, contrairement à la distance euclidienne ne mesure pas la distance géométrique entre les points en ligne droite, mais via un trajet composé de segments de droites parallèles aux axes. C’est comme si la distance euclidenne reliait les points à vol d’oiseau, alors qu’avec la distance de Manhattan, nous devions contourner les blocs de maisons du quartier pour aller d’un point A à un point B (d’où le nom de cette métrique). Elle se calcule comme suit :

\[\mathrm{D_{Manhattan}}_{j,k}=\sum_{i=1}^{n}|y_{ij}-y_{ik}|\]

Dans SciViews R, nous utiliserons dissimilarity(DF, method = "manhattan"). Ici aussi, seul l’indice de dissimilarité est défini. L’indice de similarité complémentaire n’existe pas car la valeur de l’indice de dissimlarité n’est pas borné à droite et peut varier de zéro (dissimilarité nulle, cela signifie alors que les deux individus sont identiques) à l’infini pour une différence maximale.

5.2.3 Utilisation des indices

• Les distances euclidienne ou de Manhattan sont à préférer pour les mesures environnementales ou de manière générale pour les variables quantitatives continues.

• Les distances de Bray-Curtis ou Canberra sont meilleure pour les dénombrements d’espèces (ou à chaque fois que les double zéros ne portent aucune information utile). Il s’agit souvent de variables quantitatives discrètes prenant des valeurs nulles ou positives.

À vous de jouer !

5.2.4 Propriétés des indices

Les indices varient en 0 et 1 (0 et 100%), mais les distances sont utilisées aussi comme indices de dissimilarité et peuvent varier entre 0 et \(+\infty\).

Un indice est dit métrique si il az les propriétés suivantes :

• Minimum 0 : \(I_{j, k} = 0\) si \(j = k\)

• Positif : \(I_{j, k}>0\) si \(j \neq k\)

• Symétrique : \(I_{j, k}=I_{k, j}\)

• Inégalité triangulaire : \(I_{j, k} + I_{k, l} >= I_{j, l}\)

La dernière propriété d’inégalité triangulaire est la plus difficile à obtenir, et n’est pas toujours nécessaire. Nous pouvons montrer que certains indices qui ne respectent pas cette dernière propriété sont pourtant utiles dans le contexte. Nous dirons alors d’un indice que c’est une semi-métrique s’il répond à toutes les conditions sauf la quatrième. Enfin, un indice est dit non métrique dans tous les autres cas. Le tableau suivant reprend les métriques que nous avons vues jusqu’ici, et rajoute d’autres candidats potentiels (la distance Chi carré, l’indice de corrélation ou de variance/covariance) en indiquant leur type :

Distance	Type
Bray-Curtis	semi-métrique
Canberra	métrique
Euclidienne	métrique
Manhattan	métrique
Chi carré	métrique
(corrélation)	(non métrique)
(variance/covariance)	(non métrique)

Pour en savoir plus

• Vous pouvez aussi transposer le tableau pour calculer la distance entre ses colonnes en utilisant l’argument transpose = TRUE dans dissimilarity(). Par exemple, dans le cas d’un tableau “espèces - station” (dénombrement d’espèces en différentes stations), nous pouvons comparer les stations du point de vue de la composition en espèces, mais nous pouvons aussi comparer les espèces du point de vue de leur répartition entre les stations. Pour passer d’un calcul à l’autre, nous transposerons donc le tableau (les colonnes deviennent les lignes et inversement).

• Pour bien comprendre la logique sous-jacente aux indices, il est utile de comprendre leurs équations. Si elles sont pour vous trop abstraites, une façon efficace de comprendre consiste à faire le calcul à la main. Par exemple dans le cas de l’indice de Canberra, la notion de nombre de données non nulles \(nz\) n’est pas évident. Effectuons un calcul à la main détaillé sur le tableau fictif suivant concernant trois espèces A, B, et C dénombrées en trois stations sta1, sta2 et sta3 dans le tableau nommé ex1 :

  	A	B	C
  sta1	4	0	2
  sta2	3	0	10
  sta3	1	8	0

Pour rappel, la dissimilarité de Canberra se calcule comme suit :

\[\mathrm{D_{Canberra}}_{j,k}=\frac{1}{nz}\sum_{i'=1}^{nz}\frac{\left|y_{i'j}-y_{i'k}\right|}{y_{i'j}+y_{i'k}}\]

où :

• nz est le nombre d’observations non nulles simultanément dans les deux vecteurs comparés (les doubles zéros ne sont pas pris en compte)
• i’ est l’itérateur sur toutes les valeurs non double zéros

Voici le détail du calcul (notez bien comment le double zéro pour l’espèce B entre les stations sta1 et sta2 est pris en compte dans le calcul) :

sta1_2 <- (1/2) * ((abs(4 - 3)) / (4 + 3) + (abs(2 - 10)) / (2 + 10))
round(sta1_2, 2)

# [1] 0.4

sta1_3 <- (1/3) * (abs(4 - 1) / (4 + 1) + abs(0 - 8) / (0 + 8) +
  abs(2 - 0) / (2 + 0))
round(sta1_3, 2)

# [1] 0.87

sta2_3 <- (1/3) * (abs(3 - 1) / (3 + 1) + abs(0 - 8) / (0 + 8) +
  abs(10 - 0) / (10 + 0))
round(sta2_3, 2)

# [1] 0.83

La matrice finale est la suivante :

	sta1	sta2
sta2	0.4
sta3	0.87	0.83

Vérifions en laissant R faire le calcul :

(dissimilarity(ex1, method = "canberra"))

# Dissimilarities :
#           sta1      sta2
# sta2 0.4047619          
# sta3 0.8666667 0.8333333
# 
# Metric :  
# Number of objects : 3

---

12. Le nombre de paires uniques et distinctes (pas j, j ou k, k) possibles parmi n items est \(n(n-1)/2\), soit ici pour 1262 éléments nous avons 795.691 paires.↩︎