Introduction

These are codes I used to run NMDS on IBI and metric data calculated using the MBSStools package.

Load the package

Package vegan is the package in R that you can use to run multivariate analyses (Oksanen et al., 2022).

library(vegan)

Load data

Data can be found here.

library(readxl)
metric_ibi_allsites <- read_excel("F:/GitHub Projects/thesis_codes/Part b (09-19)/NMDS (Metrics and IBI) - all sites/metrics_ibi_allsites.xlsx")

Make a matrix from IBI and metric data

ibi_wo_su <- metric_ibi_allsites[,5:ncol(metric_ibi_allsites)] #exclude the first few columns that contains SU information

Run NMDS

set.seed(1)
nmds_ibi <- metaMDS(ibi_wo_su,distance="bray",maxit=999,trymax = 500,wascores = T,k=2,autotransform = F) #gotta put autotransform on FALSE so R won't transform your data when you don't want it to.

Add back SU columns

library(dplyr) #load the package
data.scores.site <- as.data.frame(scores(nmds_ibi)$sites) #extract NMDS scores
data.scores.site.dplyr <- data.scores.site %>% mutate(Site = metric_ibi_allsites$Site,Year=metric_ibi_allsites$Year,
                                            Subwatershed=metric_ibi_allsites$Subwatershed)

Plot the graph using ggplot2

library(ggplot2) #load the package

data.scores.site.dplyr$Subwatershed = factor(data.scores.site.dplyr$Subwatershed,levels=c("B9 (2%)","B7 (12%)",
                                                                                "B9A (16%)","B1 (21%)"))

ibi_part2_allsites <- ggplot(data.scores.site.dplyr, aes(x = NMDS1, y = NMDS2)) + 
    geom_point(size = 4, stroke=2,aes( shape = as.factor(Year), fill = Subwatershed,color=Subwatershed))+
    scale_shape_manual(values = c(0,1,2,5,6,7,21,22,23,24,25)) +
  # geom_polygon(data=data.scores.dplyr,aes(x=NMDS1,y=NMDS2,group=Subwatershed),alpha=0.30); this is to add a little convex polygon to visualize the clusters better. You can try to see how it looks. 
    theme(axis.text.y = element_text(colour = "black", size = 25, face = "bold"), 
    axis.text.x = element_text(colour = "black", face = "bold", size = 25), 
    legend.text = element_text(size = 22, face ="bold", colour ="black"), 
    legend.position = "right", axis.title.y = element_text(face = "bold", size = 18), 
    axis.title.x = element_text(face = "bold", size = 25, colour = "black"), 
    axis.title.y.left = element_text(face = "bold", size = 25, colour = "black"),
    legend.title = element_text(size = 22, colour = "black", face = "bold"), 
    panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
    legend.key=element_blank(),
    plot.title = element_text(color = "black", size = 30, face = "bold", hjust = 0.5)) + 
    labs(
    title = "NMDS graph of Part B IBI and Metrics") + 
    theme(axis.title.x = element_text(margin=margin(t=10)), #add margin to x-axis title
        axis.title.y = element_text(margin=margin(r=10)))+
    labs(x = "NMDS1", fill = "Sub-watershed",color="Sub-watershed", y = "NMDS2", shape = "Year")  + 
    scale_fill_manual(values = c( "#ffeea1","#feb14c","#fc4d2a","#b10026")) +
    scale_color_manual(values = c( "#ffeea1","#feb14c","#fc4d2a","#b10026")) +
    geom_text(aes(label=Site),hjust=0.55, vjust=1.5,size=7) + 
    geom_hline(yintercept=0, linetype="dashed", 
                color = "black", size=1) + #add horizontal and vertical lines at 0
    geom_vline(xintercept=0, linetype="dashed", 
                color = "black", size=1)

ibi_part2_allsites

ggsave("Part 2 IBI NMDS - all sites.jpeg",ibi_part2_allsites,width=25,height=13)

Plot a biplot

The job now is to add some arrows that correspond to the “metric (or species)” scores extract from above.

data.scores.metric <- as.data.frame(scores(nmds_ibi)$species) #extract NMDS metric scores

ibi_allsite_biplot <- ibi_part2_allsites + 
  geom_segment(data=data.scores.metric, aes(x=0, xend=NMDS1, y=0, yend=NMDS2), 
               color="black", arrow=arrow(length=unit(0.01,"npc"))) + #add arrow
  geom_text(data=data.scores.metric, 
            aes(x=NMDS1,y=NMDS2,label=rownames(data.scores.metric),
                hjust=0.35*(1.5-sign(NMDS1)),vjust=0.75*(1-sign(NMDS2))), 
            color="black", size=6) + #add text to arrow
  geom_hline(yintercept=0, linetype="dashed", 
                color = "black", size=1) + #add horizontal and vertical lines at 0
  geom_vline(xintercept=0, linetype="dashed", 
                color = "black", size=1)

ibi_allsite_biplot

ggsave("NMDS Metric All Sites Biplot.jpeg", ibi_allsite_biplot,width=25,height=12)

Plot a Screeplot

library(goeveg) #load package goeveg

dimcheckMDS(ibi_wo_su,distance = "bray",k=10,trymax=500,autotransform = F)

Perform a Mantel test

The purpose of this step is to examine the correlation between between the two distance (or dissimilarity) matrices. Biologically speaking, I wanted to investigate whether or not doing NMDS on count data is different from doing so on IBI. Thus, it might be valuable later on to choose which type of data to go with to get the most information.

Load data

Data can be found here.

library(readxl)
count_combined_allsites <- read_excel("F:/GitHub Projects/thesis_codes/Part b (09-19)/NMDS (Count) - all sites - spatial question/count_combined_allsites.xlsx", 
    sheet = "CHIRAE F OLIGTA O")

Clean up data

Since the data set contains a lot of columns with zeros, it is essential to clean up all zeros to have a clean data set.

count_nozero <- count_combined_allsites[, colSums(count_combined_allsites != 0) > 0]

#Another way is to use dplyr 

library(dplyr)
count_nozero_dplyr <- count_combined_allsites %>% select(where(~ any(. != 0))) #this will yield the same file

#Write out a non-zero file to use for later. I will use PC-ORD v9 to cross check NMDS outputs
library("writexl")

write_xlsx(count_nozero,"F:/One Drive - Towson U/Thesis stuff/Data Analysis/Part b (2009 - recent)/count_nozero_allsites.xlsx")

Make a matrix from count data

library(funrar)
count_wo_su <- count_nozero[,6:ncol(count_nozero)] #exclude the first few columns that contains SU information
rel_count_wo_su <- make_relative(as.matrix(count_nozero[,6:ncol(count_nozero)]))
library(vegan)

###Calculate dissimilarity matrices

diss_ibi <- vegdist(ibi_wo_su, method = "bray")
diss_count <- vegdist(count_wo_su, method = "bray") 
mantel(diss_count,diss_ibi,method="pearson", permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = diss_count, ydis = diss_ibi, method = "pearson",      permutations = 10000) 
## 
## Mantel statistic r: 0.5792 
##       Significance: 9.999e-05 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0859 0.1101 0.1346 0.1633 
## Permutation: free
## Number of permutations: 10000