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_combined <- read_excel("F:/GitHub Projects/thesis_codes/Part b (09-19)/NMDS (Metrics and IBI) - B1, 7, 10, REF 12-19/metric_ibi_combined.xlsx", 
    sheet = "Metric & IBI")

Make a matrix from IBI and metric data

ibi_wo_su <- metric_ibi_combined[,4:ncol(metric_ibi_combined)] #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_combined$Site,Year=metric_ibi_combined$Year,
                                            Subwatershed=metric_ibi_combined$Subwatershed)

Plot the graph using ggplot2

library(ggplot2) #load the package

ibi_part2_graph <- ggplot(data.scores.site.dplyr, aes(x = NMDS1, y = NMDS2)) + 
    geom_point(size = 4,aes( shape = as.factor(Year), fill = Subwatershed,color=Subwatershed))+
    scale_shape_manual(values = c(11,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("#b10026", "#feb14c","#ffeea1","#fc4d2a")) +
    scale_color_manual(values = c("#b10026", "#feb14c","#ffeea1","#fc4d2a")) +
    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_graph

ggsave("Part 2 IBI NMDS (B1, 7, 10, REF).jpeg",ibi_part2_graph,width=25,height=12)

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

nmds_metric_part2_biplot <- ibi_part2_graph + 
  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.55*(1-sign(NMDS1)),vjust=0.5*(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)

nmds_metric_part2_biplot

ggsave("NMDS Metric Part B Biplot.jpeg", nmds_metric_part2_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)