Science des données biologiques 2

Réalisé par le service d'Écologie numérique des Milieux aquatiques, Université de Mons (Belgique)

Préambule

Si vous n’avez jamais utilisé de tutoriel “learnr”, familiarisez-vous d’abord avec son interface ici.

Ne vous trompez pas dans votre adresse mail et votre identifiant Github

N’oubliez pas de soumettre votre réponse après chaque exercice

Conformément au RGPD (Règlement Général sur la Protection des Données), nous sommes tenus de vous informer de ce que vos résultats seront collecté afin de suivre votre progression. Les données seront enregistrées au nom de l’utilisateur apparaissant en haut de cette page. Corrigez si nécessaire ! En utilisant ce tutoriel, vous marquez expressément votre accord pour que ces données puissent être collectées par vos enseignants et utilisées pour vous aider et vous évaluer. Après avoir été anonymisées, ces données pourront également servir à des études globales dans un cadre scientifique et/ou éducatif uniquement.

Contexte

Des scientifiques ont réalisés des mesures sur l’Île Barro Colorado, Panama. Il s’agit d’une ile artificielle située sur le lac Gatùn.

Le jeu de données se nomme BCI et provient du package vegan. Pour plus de simplicité nous le nommerons bci.

bci <- read("BCI", package = "vegan")

MDS simplifiée sous SciViews::R

Afin de vous simplifier la réalisation des mds, des fonctions sont mises à votre disposition.

SciViews::R()
library(broom)

# function mds for several multidimensionnal scaling functions ------
mds <- function(dist, k = 2, type = c("metric", "nonmetric", "cmdscale",
"wcmdscale", "sammon", "isoMDS", "monoMDS", "metaMDS"), p = 2, ...) {
  type <- match.arg(type)
  res <- switch(type,
    metric = ,
    wcmdscale = structure(vegan::wcmdscale(d = dist, k = k, eig = TRUE, ...),
      class = c("wcmdscale", "mds", "list")),
    cmdscale = structure(stats::cmdscale(d = dist, k = k, eig = TRUE, ...),
      class = c("cmdscale", "mds", "list")),
    nonmetric = ,
    metaMDS = structure(vegan::metaMDS(comm = dist, k = k, ...),
      class = c("metaMDS", "monoMDS", "mds", "list")),
    isoMDS = structure(MASS::isoMDS(d = dist, k = k, ...),
      class = c("isoMDS", "mds", "list")),
    monoMDS = structure(vegan::monoMDS(dist = dist, k = k, ...),
      class = c("monoMDS", "mds", "list")),
    sammon = structure(MASS::sammon(d = dist, k = k, ...),
      class = c("sammon", "mds", "list")),
    stop("Unknown 'mds' type ", type)
  )
  # For non-metric MDS, we add also data required for the Shepard plot
  if (type %in% c("nonmetric", "sammon", "isoMDS", "monoMDS", "metaMDS"))
    res$Shepard <- MASS::Shepard(d = dist, x = res$points, p = p)
    res
}
class(mds) <- c("function", "subsettable_type")

# plot.mds : MDS2 ~ MDS1 --------------------------------
plot.mds <- function(x, y, ...) {
  points <- tibble::as_tibble(x$points, .name_repair = "minimal")
  colnames(points) <- paste0("mds", 1:ncol(points))
  
  plot(data = points, mds2 ~ mds1,...)
}

autoplot.mds <- function(object,  labels, ...) {
  points <- tibble::as_tibble(object$points, .name_repair = "minimal")
  colnames(points) <- paste0("mds", 1:ncol(points))
  
  if (!missing(labels)) {
    if (length(labels) != nrow(points))
      stop("You must provide a character vector of length ", nrow(points),
        " for 'labels'")
    points$.labels <- labels
    chart::chart(points, mds2 ~ mds1 %label=% .labels, ...) +
      geom_point() +
      ggrepel::geom_text_repel() +
      coord_fixed(ratio = 1)
  } else {# Plot without labels
    chart::chart(points, mds2 ~ mds1, ...) +
      geom_point() +
      coord_fixed(ratio = 1)
  }
}

shepard <- function(dist, mds, p = 2)
  structure(MASS::Shepard(d = dist, x = mds$points, p = p),
    class = c("shepard", "list"))

plot.shepard <- function(x, y, l.col = "red", l.lwd = 1,
xlab = "Observed Dissimilarity", ylab = "Ordination Distance", ...) {
  she <- tibble::as_tibble(x, .name_repair = "minimal")
  
  plot(data = she, y ~ x, xlab = xlab, ylab = ylab, ...)
  lines(data = she, yf ~ x, type = "S", col = l.col, lwd = l.lwd)
}

autoplot.shepard <- function(object,  alpha = 0.5, l.col = "red", l.lwd = 1,
xlab = "Observed Dissimilarity", ylab = "Ordination Distance", ...) {
  she <- tibble::as_tibble(object)
  
  chart(data = she, y ~ x) +
    geom_point(alpha = alpha) +
    geom_step(chart::f_aes(yf ~ x), direction = "vh", col = l.col, lwd = l.lwd) +
    labs(x = xlab, y = ylab)
}

# augment.mds -------------------------------------------
augment.mds <- function(x, data, ...){
  points <- as_tibble(x$points)
  colnames(points) <- paste0(".mds", 1:ncol(points))
  bind_cols(data, points)
}

# glance.mds -------------------------------------------
glance.mds <- function(x, ...){
  if ("GOF" %in% names(x)) {# Probably cmdscale() or wcmdscale() => metric MDS
    tibble::tibble(GOF1 = x$GOF[1], GOF2 = x$GOF[2])
  } else {# Non metric MDS
    # Calculate linear and non linear R^2 from the Shepard (stress) plot
    tibble::tibble(
      linear_R2 = cor(x$Shepard$y, x$Shepard$yf)^2,
      nonmetric_R2 = 1 - sum(vegan::goodness(x)^2)
    )
  }
}

Analyse en coordonnées principales (ou MDS métrique)

Réalisez une PCoA sur le jeu de données bci avec une matrice de distance de Canberra.

# function de base
dist. <- vegan::vegdist(x, method = "bray")
mds. <- mds$metric(dist.)

autoplot(mds.)
glance(mds.)
summary(bci)
dist. <- vegan::vegdist(x, method = "bray")
mds. <- mds$metric(dist.)
autplot(mds.)
glance(mds.)
dist. <- vegan::vegdist(bci, method = "canberra")
mds. <- mds$metric(dist.)
autoplot(mds.)
glance(mds.)
# TODO

MDS non métrique

# function de base
dist. <- vegan::vegdist(x, method = "bray")
mds. <- mds$nonmetric(dist.)

autoplot(mds.)
glance(mds.)
summary(bci)
dist. <- vegan::vegdist(x, method = "bray")
mds. <- mds$nonmetric(dist.)
autplot(mds.)
glance(mds.)
dist. <- vegan::vegdist(bci, method = "canberra")
mds. <- mds$nonmetric(dist.)
autoplot(mds.)
glance(mds.)
# TODO

Conclusion

Vous venez de terminer votre séance d’exercice.

Laissez nous vos impressions sur cet outil pédagogique ou expérimentez encore dans la zone ci-dessous. Rappelez-vous que pour placer un commentaire dans une zone de code R, vous devez utilisez un dièse (#) devant vos phrases.

# Ajout de commentaires 
# ...
# Not yet...

Guyliann Engels & Philippe Grosjean