NMDS

Script NMDS plots

Packages

library(reshape2)

library(ggplot2) 

library(dplyr) 

library(readxl) 

library(vegan) 

Data

il <- read_xlsx("Inv_Long.xlsx")

desc <- read_xlsx("Replicate_Descriptives.xlsx")

iw <- merge(desc, il, by= "Replicate_ID")

iw <- iw %>% distinct(.keep_all = TRUE)

iw <- iw %>% filter(Sample_ID != "KO_1_S_D" & Name != "Chironomidae.spp.")

1. NMDS - Abundance per Replicate ID

iwr<- dcast(iw, Replicate_ID + Code ~ Name, value.var = "Count",
            fun.aggregate = sum,fill = 0)

Determine number of dimensions to use in the NMDS

Note that more dimensions = higher complexity to interpret.

Create distance matrix.

dist_rep <- vegdist(iwr[, 3:35], method = "bray")  

Function that performs 1-10 dimensions NMDS & plot stress levels.

NMDS.scree <- function(x) { plot(rep(1, 10), replicate(10, 
          metaMDS(x, autotransform = F, k = 1)$stress),  
          xlim = c(1, 10),ylim = c(0, 0.30), xlab = "# of Dimensions",  
          ylab = "Stress", main = "NMDS stress plot") 
  for (i in 1:10) {  points(rep(i + 1,10),replicate(10, 
          metaMDS(x, autotransform = F, k = i + 1)$stress)) } }

Plot NMDS stress for my distance matrix.

Chosen NMDS

nmds_iwr <- metaMDS(comm = iwr[ , 3:35], k = 3,distance = "bray", trymax = 100)

print(nmds_iwr)

## 
## Call:
## metaMDS(comm = iwr[, 3:35], distance = "bray", k = 3, trymax = 100) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(iwr[, 3:35])) 
## Distance: bray 
## 
## Dimensions: 3 
## Stress:     0.1545876 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 37 tries
## The best solution was from try 18 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(iwr[, 3:35]))'

Lowest stress: 0.15, type 1: (can be useful but has potential to mislead)

Evaluate Fit

stressplot(nmds_iwr)

Good Non-Metric fit, R2= 0.976

Extract NMDS scores

data.iwr = as.data.frame(scores(nmds_iwr)$sites)

species.iwr = as.data.frame(scores(nmds_iwr)$species)

species.iwr$species <- rownames(species.iwr)

Add descriptive variables

data.iwr$Sample = iwr$Sample_ID

data.iwr$Code = iwr$Code

Plot NMDS

ggplot() + 
  geom_text(data = species.iwr, aes(x = NMDS1, y = NMDS2, label = species), 
            alpha = 0.9, size = 3) + 
  geom_point(data = data.iwr, aes(x = NMDS1, y = NMDS2, color = data.iwr$Code),
             size = 3) +  
  annotate(geom = "label", x = 0.5, y = 1.25, size = 5,fill="lightgrey",
           label = paste("Stress: ", round(nmds_iwr$stress, digits = 3))) +
  theme_bw() + 
  theme(legend.position = "right", text = element_text(size = 15) ) + 
  labs(colour="Site")

2. NMDS - Mean abundance per Sample ID —-

Calculate mean abundance per by Sample_ID.

iwm <- dcast(iw, Sample_ID + Code ~ Name, value.var = "Count", 
             fun.aggregate = mean,fill = 0)

Determine number of dimensions to use in the NMDS

dist_sam <- vegdist(iwm[, 3:35], method = "bray") 

Plot NMDS stress for my distance matrix.

Chosen NMDS

nmds_iwm <- metaMDS(comm = iwm[ , 3:35], k = 3,distance = "bray", trymax = 100)

print(nmds_iwm) 

## 
## Call:
## metaMDS(comm = iwm[, 3:35], distance = "bray", k = 3, trymax = 100) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(iwm[, 3:35])) 
## Distance: bray 
## 
## Dimensions: 3 
## Stress:     0.1301082 
## Stress type 1, weak ties
## Best solution was repeated 2 times in 20 tries
## The best solution was from try 16 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(iwm[, 3:35]))'

Evaluate Fit

stressplot(nmds_iwm) 

Good Non-Metric fit, R2= 0.983

Extract NMDS scores

data.iwm = as.data.frame(scores(nmds_iwm)$sites)

species.iwm = as.data.frame(scores(nmds_iwm)$species)

species.iwm$species <- rownames(species.iwm)

Add descriptive variables

data.iwm$Sample = iwm$Sample_ID

data.iwm$Code = iwm$Code

Plot NMDS

ggplot() + 
  geom_text(data = species.iwm, aes(x = NMDS1, y = NMDS2, label = species), 
            alpha = 0.9, size = 3) + 
  geom_point(data = data.iwm, aes(x = NMDS1, y = NMDS2, color = data.iwm$Code),
             size = 3) +  
  annotate(geom = "label", x = 0.5, y = 1.25, size = 5,fill="lightgrey", 
           label = paste("Stress: ", round(nmds_iwm$stress, digits = 3))) +
  theme_bw() + 
  theme(legend.position = "right", text = element_text(size = 15) ) +
  labs(colour="Site")