non-Metric Multidimentional Scaling (nMDS)

Step 1

Load libraries.
Cargar las librerias que necesitas.

libraries <- c("vegan", "ggplot2", "dplyr")
lapply(libraries, require, character.only = TRUE)

## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE

Step 2

Data
Cargar los datos.

laselva_full  <- read.csv("D:/Curriculum/07_ Cursos/Course_Multivariate_Stat_for_Ecological_Data/data/nMDS_laSelva.csv")
head(laselva_full)

Step 3

Seleccionar la matrix de especies.

laselva_sp <- dplyr::select(laselva_full, Baetid:Collemb)
str(laselva_sp)
ncol(laselva_sp)
nrow(laselva_sp)

non-Metric Multidimentional Scaling (nMDS)

Step 1: nMDS

Ejecutar un nMDS

set.seed(14) # Con este comando, siempre comenzara del mismo lugar.
laselva.mds <- metaMDS(laselva_sp, distance = "bray", k = 2, trymax=100)  #using all the defaults

Step 2: nMDS

Resumen de los reusltados

laselva.mds

## 
## Call:
## metaMDS(comm = laselva_sp, distance = "bray", k = 2, trymax = 100) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(laselva_sp)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2139587 
## Stress type 1, weak ties
## No convergent solutions - best solution after 100 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(laselva_sp))'

Lo más importante de ver en este “output” es el stress

As a rule of thumb literature has identified the following cut-ff values for stress-level:
> 0.2 is poor (risks for false interpretation).
0.1 - 0.2 is fair (some distances can be misleading for interpretation).
0.05 - 0.1 is good (can be confident in inferences from plot).
< 0.05 is excellent (this can be rare).

Step 3: nMDS

Plot the results.

plot(laselva.mds, type="t")

Step 4: nMDS

Stressplot. Large scatter around the line suggests that original dissimilarities are not well preserved in the reduced number of dimensions. Looks pretty good in this case.

stressplot(laselva.mds)

Step 5: nMDS

Otra metrica para evaluar el ajuste (buen desempeño) del nMDS es mirar el “Goodness of fit”.

gof <- goodness(laselva.mds)
# gof

{plot (laselva.mds, display = 'sites', type = 't', main = 'Goodness of fit') # this function draws NMDS ordination diagram with sites
points (laselva.mds, display = 'sites', cex = 2*gof/mean(gof))} # and this adds the points with size reflecting goodness of fit (bigger = worse fit)

Step 6: nMDS

Calcular el % de variacion explicado por cada eje

MDS <- cmdscale(vegdist(laselva_sp, method = "bray"), k = 2, eig = T, add = T )
round(MDS$eig*100/sum(MDS$eig),1)

##   [1] 25.6 10.5  4.8  3.9  3.4  3.2  2.8  2.6  2.4  2.1  1.9  1.7  1.7  1.4  1.2
##  [16]  1.2  1.1  1.0  1.0  1.0  0.9  0.9  0.8  0.8  0.8  0.7  0.7  0.7  0.6  0.6
##  [31]  0.6  0.5  0.5  0.5  0.5  0.5  0.4  0.4  0.4  0.4  0.4  0.4  0.4  0.4  0.3
##  [46]  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.2  0.2
##  [61]  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2
##  [76]  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.2  0.1  0.1  0.1  0.1
##  [91]  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1  0.1
## [106]  0.0  0.0  0.0  0.0

Visualizar los datos con `ggplot2`

Step 1: Data visualization

Score son los axis de la ordenacion

data.scores <- as.data.frame(scores(laselva.mds, "sites")) 
head(data.scores)

##        NMDS1       NMDS2
## 1 -0.3467815  0.01961633
## 2 -0.1380800 -0.07419054
## 3 -0.3924111 -0.07021848
## 4  0.0671686 -0.12208157
## 5 -0.5214802  0.02659047
## 6 -0.5771004  0.24303366

Step 2: Data visualization

Agregar las especies al dataframe llamado data.scores

species.scores <- as.data.frame(scores(laselva.mds, "species"))  #Using the scores function from vegan to extract the species scores and convert to a data.frame
#species.scores$species <- rownames(species.scores)  # create a column of species, from the rownames of species.scores
head(species.scores)  #look at the data

##              NMDS1       NMDS2
## Baetid  -0.2396226  0.34625957
## America -0.9641251 -0.03831084
## Calliba -0.1270288  0.05090964
## Guajir  -1.3459779 -0.11980049
## Caenis  -0.3220985 -0.05372070
## Ulmerit -1.1376205 -0.02146410

Step 3: Data visualization

Extraer los meses de la base de datos original (laselva_full) y agregarla a data.scores

month  <- dplyr::select(laselva_full, month)
data.scores$month <- unlist(month)
head(data.scores)

##        NMDS1       NMDS2 month
## 1 -0.3467815  0.01961633   one
## 2 -0.1380800 -0.07419054   one
## 3 -0.3924111 -0.07021848   one
## 4  0.0671686 -0.12208157   one
## 5 -0.5214802  0.02659047   two
## 6 -0.5771004  0.24303366   two

Step 4: Data visualization

Hay cuatro grupos, previamente identificados por un cluster

grp.a <- data.scores[data.scores$month == "one", ][chull(data.scores[data.scores$month == 
          "one", c("NMDS1", "NMDS2")]),]  # hull values for grp A
grp.b <- data.scores[data.scores$month == "two", ][chull(data.scores[data.scores$month == 
          "two", c("NMDS1", "NMDS2")]),]  # hull values for grp A  
grp.c <- data.scores[data.scores$month == "three", ][chull(data.scores[data.scores$month == 
          "three", c("NMDS1", "NMDS2")]),]  # hull values for grp A  
grp.d <- data.scores[data.scores$month == "four", ][chull(data.scores[data.scores$month == 
          "four", c("NMDS1", "NMDS2")]),]  # hull values for grp A  

hull.data <- rbind(grp.a, grp.b, grp.c, grp.d)
head(hull.data)

##         NMDS1       NMDS2 month
## 4   0.0671686 -0.12208157   one
## 3  -0.3924111 -0.07021848   one
## 1  -0.3467815  0.01961633   one
## 35  0.4702855 -0.07926859   two
## 37  0.4610409 -0.17262389   two
## 32  0.3635824 -0.22033234   two

Step 5: Data visualization

ggplot() + 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=month,group=month),alpha=0.30) + # add the convex hulls
#  geom_text(data=species.scores,aes(x=NMDS1,y=NMDS2,label=species),alpha=0.5) +  # add the species labels
  geom_point(data=data.scores,aes(x=NMDS1,y=NMDS2,shape=month,colour=month),size=4) + # add the point markers
#  scale_colour_manual(values=c("one" = "red", "two" = "blue",
#                                "three"="black", "four"="green")) +
  coord_equal() +
  theme_bw() + 
  theme(axis.text.x = element_blank(),  # remove x-axis text
        axis.text.y = element_blank(), # remove y-axis text
        axis.ticks = element_blank(),  # remove axis ticks
        axis.title.x = element_text(size=18), # remove x-axis labels
        axis.title.y = element_text(size=18), # remove y-axis labels
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())

non-Metric Multidimentional Scaling (nMDS)

Pablo E. Gutiérrez-Fonseca

29 September, 2022

Step 1

Step 2

Step 3

non-Metric Multidimentional Scaling (nMDS)

Step 1: nMDS

Step 2: nMDS

Step 3: nMDS

Step 4: nMDS

Step 5: nMDS

Step 6: nMDS

Visualizar los datos con `ggplot2`

Step 1: Data visualization

Step 2: Data visualization

Step 3: Data visualization

Step 4: Data visualization

Step 5: Data visualization

non-Metric Multidimentional Scaling (nMDS)

Pablo E. Gutiérrez-Fonseca

29 September, 2022

Step 1

Step 2

Step 3

non-Metric Multidimentional Scaling (nMDS)

Step 1: nMDS

Step 2: nMDS

Step 3: nMDS

Step 4: nMDS

Step 5: nMDS

Step 6: nMDS

Visualizar los datos con ggplot2

Step 1: Data visualization

Step 2: Data visualization

Step 3: Data visualization

Step 4: Data visualization

Step 5: Data visualization

Visualizar los datos con `ggplot2`