new learnR 06b and 06c

Author: GuyliannEngels

inst/tutorials/06b_mds/images/BioDataScience-128.png

@@ -0,0 +1,398 @@
+---
+title: "Positionnement multidimensionnel (MDS)"
+author: "Guyliann Engels & Philippe Grosjean"
+output:
+  learnr::tutorial
+tutorial:
+  id: "sdd2.06b"
+  version: 1.0.0
+runtime: shiny_prerendered
+---
+
+```{r setup, include=FALSE}
+# Packages ----
+library(learnr)
+library(knitr)
+SciViews::R()
+library(BioDataScience)
+library(ade4)
+
+# Collect informations ------
+options(tutorial.event_recorder = BioDataScience::record_sdd)
+tutorial_options(exercise.checker = BioDataScience::checker_sdd)
+tutorial_options(exercise.timelimit = 60)
+tutorial_options(exercise.cap = "Code R")
+knitr::opts_chunk$set(echo = FALSE, comment = NA)
+
+# Preparation dataset ------
+bci <- read("BCI", package = "vegan")
+# sciviews functions ---------------
+
+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)
+    )
+  }
+}
+```
+
+```{r, echo=FALSE}
+fixedRow(
+  column(9, div(
+    img(src = 'images/BioDataScience-128.png', align = "left"),
+    h1("Science des données biologiques 2"),
+    "Réalisé par le service d'Écologie numérique des Milieux aquatiques, Université de Mons (Belgique)"
+  )),
+  column(3, div(
+    textInput("user", "Utilisateur :", ""),
+    textInput("email", "Email :", "")
+  ))
+)
+textOutput("user") # This is newer shown, but required to trigger an event!
+textOutput("email") # Idem!
+```
+
+```{r, context="server"}
+output$user <- renderText({BioDataScience::user_name(input$user);""})
+output$email <- renderText({BioDataScience::user_email(input$email);""})
+updateTextInput(session, "user", value = BioDataScience::user_name())
+updateTextInput(session, "email", value = BioDataScience::user_email())
+```
+
+## Préambule
+
+Si vous n'avez jamais utilisé de tutoriel "learnr", familiarisez-vous d'abord avec son interface [ici](http://biodatascience-course.sciviews.org/sdd-umons/learnr.html).
+
+![](images/attention.jpg)
+
+**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](https://ec.europa.eu/info/law/law-topic/data-protection/reform/rules-business-and-organisations/principles-gdpr_fr)), 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`. 
+
+```{r, echo = TRUE}
+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. 
+
+```{r, echo = TRUE}
+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`. 
+
+```{r, echo = TRUE, eval = FALSE}
+# function de base
+dist. <- vegan::vegdist(x, method = "bray")
+mds. <- mds$metric(dist.)
+
+autoplot(mds.)
+glance(mds.)
+```
+
+```{r mds, exercise = TRUE}
+summary(bci)
+```
+
+```{r mds-hint-1}
+dist. <- vegan::vegdist(x, method = "bray")
+mds. <- mds$metric(dist.)
+autplot(mds.)
+glance(mds.)
+```
+
+```{r mds-hint-2}
+dist. <- vegan::vegdist(bci, method = "canberra")
+mds. <- mds$metric(dist.)
+autoplot(mds.)
+glance(mds.)
+```
+
+```{r mds-check}
+# TODO
+```
+
+```{r qu_mds}
+question("Considérez vous que cette PCoA est une bonne PCoA ?",
+           answer("oui, sur base des valeurs de goodness-of-fit, 
+                   nous pouvons observer que le PCoA exprime une 
+                   grande part de la variance"),
+           answer("non, sur base des valeurs de goodness-of-fit, 
+                   nous pouvons observer que le PCoA n'exprime pas 
+                   une grande part de la variance. Il est préférable 
+                   de réaliser une MDS non métrique.", 
+                  correct = TRUE),
+          allow_retry = TRUE)
+```
+
+## MDS non métrique
+
+```{r, echo = TRUE, eval = FALSE}
+# function de base
+dist. <- vegan::vegdist(x, method = "bray")
+mds. <- mds$nonmetric(dist.)
+
+autoplot(mds.)
+glance(mds.)
+```
+
+```{r nmds, exercise = TRUE}
+summary(bci)
+```
+
+```{r nmds-hint-1}
+dist. <- vegan::vegdist(x, method = "bray")
+mds. <- mds$nonmetric(dist.)
+autplot(mds.)
+glance(mds.)
+```
+
+```{r nmds-hint-2}
+dist. <- vegan::vegdist(bci, method = "canberra")
+mds. <- mds$nonmetric(dist.)
+autoplot(mds.)
+glance(mds.)
+```
+
+```{r nmds-check}
+# TODO
+```
+
+```{r qu_nmds}
+question("Considérez vous que cette MDS non métrique est une 
+          MDS non métrique de qualité?",
+           answer("oui, sur base de la valeur du R^2 linéaire élevée, 
+                   nous pouvons considérer que la MDS non métrique 
+                   est exprime une grande part de la variance"),
+           answer("oui, sur base des valeurs de R^2 linéaire et de R^2
+                   non métrique élevées, nous pouvons considérer que la
+                   MDS non métrique est de qualité"),
+          answer("oui, sur base de la valeur du R^2 non métrique élevée,
+                  nous pouvons considérer que la MDS non métrique 
+                  est de qualité", correct = TRUE),
+           answer("non, sur base des valeurs R^2 linéaire et de R^2
+                   non métrique élevées, nous pouvons considérer 
+                   que la MDS non métrique n'est pas optimale."),
+          allow_retry = TRUE)
+```
+
+
+## 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.
+
+```{r comm, exercise=TRUE, exercise.lines = 8}
+# Ajout de commentaires 
+# ...
+```
+
+```{r comm-check}
+# Not yet...
+```