Electrofishing NMDS Analysis for 2022 Shiawassee Master’s Project Team

Load packages

rm(list = ls())
# Install
if(!require("vegan")) install.packages("vegan")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("ggplot2")) install.packages("ggplot2")

# Load
library(vegan)
library(tidyverse)
library(ggplot2)

Set up environment and output language

# Set your Working Directory to where your .csv file(s) are located.
# You can also do this manually from th RStudio toolbar: 
# Session > Set Working Directory > Choose Directory...
Sys.setenv(LANG = "en")

#### Changed
setwd("C:/Users/hayde/OneDrive/Desktop/UMichigan/Masters project/code and anlysis/Electrofishing/PERMANOVA/drafts")


#{r "C:/Users/hayde/OneDrive/Desktop/UMichigan/Masters project/code and anlysis/Electrofishing/PERMANOVA/final version/ef2_perm", include=FALSE}

#file2 <- read.csv('C:/Users/hayde/OneDrive/Desktop/UMichigan/Masters project/code and anlysis/Electrofishing/PERMANOVA/final version/ef2_perm.csv')

Loading data, setup

#read in your .csv files
EFNMDS = read.csv("C:/Users/hayde/OneDrive/Desktop/UMichigan/Masters project/code and anlysis/Electrofishing/PERMANOVA/drafts/EF_permanova.csv", header = T)

# site, month, veg zone, LOTU (as column headers)

# changing NA values to zeroes
EFNMDS[is.na(EFNMDS)] = 0

# changing numerical months to nominal months
EFNMDS$month[EFNMDS$month == "6"] = "June"
EFNMDS$month[EFNMDS$month == "7"] = "July"
EFNMDS$month[EFNMDS$month == "8"] = "August"

# editing veg zone names
EFNMDS$Vegetation.Type[EFNMDS$Vegetation.Type == "Open Water"] = "Open"
EFNMDS$Vegetation.Type[EFNMDS$Vegetation.Type == "Submerged Aquatic Vegetation"] = "SAV"

# subsets of SiteID, parameter of interest, LOTU counts
EFNMDS_unit = subset(EFNMDS[,c(1, 3, 5:36)])
EFNMDS_month = subset(EFNMDS[,c(1, 2, 5:36)])
EFNMDS_veg = subset(EFNMDS[,c(1, 4:36)])

Calculate the distance matrix.

# Convert cpue data from 'UMSEAS-SNWR-2022_EF-NMDS-PERMANOVA.csv' to relative abundance and Calculate user-defined distance matrix using "method=Bray-Curtis Dissimilarity".
# Note: Result matrix from this block will be later used in NMDS.
EF_rel = decostand(EFNMDS[,5:36], method = "total")
EF_distmat = vegdist(EF_rel, method = "bray")

# Read in the distance matrix from the prior function as a matrix and Write it back to your working directory.
EF_distmat = as.matrix(EF_distmat, labels = T)
write.csv(EF_distmat, "UMSEAS-SNWR-2022_EFNMDS_braydistancematrix.csv")

Perform NMDS using the full IVI Matrix and the Bray-Curtis Dissimilarity Matrix and plot.

# Note: The metaMDS() function uses random starts, but then performs PCA rotation on the results so that axis 1 (i.e., 'MDS1') represents the highest percent variation in your data; Look for stress <0.2
EF_NMDS = metaMDS(EF_distmat, distance = "bray", k = 3, maxit = 999, trymax = 500, wascores = TRUE)
## Run 0 stress 0.1267707 
## Run 1 stress 0.1276656 
## Run 2 stress 0.1254376 
## ... New best solution
## ... Procrustes: rmse 0.08538499  max resid 0.3201498 
## Run 3 stress 0.1267751 
## Run 4 stress 0.1267717 
## Run 5 stress 0.1284168 
## Run 6 stress 0.1284151 
## Run 7 stress 0.1257297 
## ... Procrustes: rmse 0.06467532  max resid 0.2041401 
## Run 8 stress 0.1284211 
## Run 9 stress 0.1243046 
## ... New best solution
## ... Procrustes: rmse 0.07455479  max resid 0.3227066 
## Run 10 stress 0.1272586 
## Run 11 stress 0.1265138 
## Run 12 stress 0.1284217 
## Run 13 stress 0.1267716 
## Run 14 stress 0.1257163 
## Run 15 stress 0.1254412 
## Run 16 stress 0.1279248 
## Run 17 stress 0.1254406 
## Run 18 stress 0.1255137 
## Run 19 stress 0.1283152 
## Run 20 stress 0.1272556 
## Run 21 stress 0.1284213 
## Run 22 stress 0.1243025 
## ... New best solution
## ... Procrustes: rmse 0.01007281  max resid 0.04326075 
## Run 23 stress 0.1267705 
## Run 24 stress 0.1264209 
## Run 25 stress 0.1257085 
## Run 26 stress 0.1257342 
## Run 27 stress 0.1242999 
## ... New best solution
## ... Procrustes: rmse 0.001045281  max resid 0.004530241 
## ... Similar to previous best
## *** Best solution repeated 1 times
goodness(EF_NMDS)
##  [1] 0.019958594 0.014080170 0.027817216 0.015525442 0.023796296 0.021804672
##  [7] 0.009582645 0.018460512 0.010950222 0.025263185 0.019938847 0.020756741
## [13] 0.023191769 0.029546762 0.017423499 0.031380505 0.022877412 0.021439411
## [19] 0.023045128 0.019061854 0.031008115 0.021619195 0.031186945 0.020301969
## [25] 0.027440808 0.023505420 0.025567774 0.020418198 0.028162638 0.014785222
EF_NMDS
## 
## Call:
## metaMDS(comm = EF_distmat, distance = "bray", k = 3, trymax = 500,      wascores = TRUE, maxit = 999) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     EF_distmat 
## Distance: user supplied 
## 
## Dimensions: 3 
## Stress:     0.1242999 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 27 tries
## The best solution was from try 27 (random start)
## Scaling: centring, PC rotation 
## Species: scores missing
stressplot(EF_NMDS)

# Create a new df that gets your MDS output data for plotting by adding the 'points' column to the EF_NMDS df.
site.scores = as_tibble(EF_NMDS$points)
site.scores
## # A tibble: 30 × 3
##        MDS1    MDS2     MDS3
##       <dbl>   <dbl>    <dbl>
##  1 -0.215   -0.141  -0.0591 
##  2 -0.141   -0.197   0.00163
##  3  0.00897 -0.150   0.122  
##  4 -0.143   -0.175  -0.0772 
##  5 -0.121   -0.0766  0.133  
##  6 -0.0956  -0.154  -0.00313
##  7 -0.0754  -0.158  -0.0309 
##  8 -0.150   -0.0859 -0.0384 
##  9 -0.0879  -0.159  -0.0784 
## 10  0.0685  -0.212  -0.0843 
## # … with 20 more rows
# Add in the grouping variables (i.e., treatments) to the 'site.scores' df to call to when color/shape coding. 
site.scores$Unit = factor(paste(EFNMDS$Unit))
site.scores$Vegetation.Type = factor(paste(EFNMDS$Vegetation.Type))
site.scores$Month = paste(EFNMDS$month)
#site.scores$Month = factor(site.scores$Month, levels= c("Jun", "Jul", "Aug"))
site.scores
## # A tibble: 30 × 6
##        MDS1    MDS2     MDS3 Unit  Vegetation.Type Month 
##       <dbl>   <dbl>    <dbl> <fct> <fct>           <chr> 
##  1 -0.215   -0.141  -0.0591  MC    Open            June  
##  2 -0.141   -0.197   0.00163 MC    Open            July  
##  3  0.00897 -0.150   0.122   MC    SAV             August
##  4 -0.143   -0.175  -0.0772  MC    SAV             June  
##  5 -0.121   -0.0766  0.133   MC    SAV             July  
##  6 -0.0956  -0.154  -0.00313 MC    SAV             August
##  7 -0.0754  -0.158  -0.0309  MC    SAV             June  
##  8 -0.150   -0.0859 -0.0384  MC    SAV             July  
##  9 -0.0879  -0.159  -0.0784  MC    SAV             August
## 10  0.0685  -0.212  -0.0843  MN    SAV             June  
## # … with 20 more rows

##Time to plot!

## MONTH
# Color code by MONTH, code shapes by Unit.

#MDS1 vs. MDS2: (The first time, run without 'coord_cartesian()', then add that customization script back in after you've seen what the base plot looks like. Coord_cartesian() makes your axes' min/max values consistent. 
# Annotate() provides a legend overlaid on the plot (at coordinates x,y) with your final k and stress values that you'll need to input from the NMDS output.
EF2022NMDS_MONTH_k3_MDS12 = ggplot(site.scores, aes(x=MDS1, y=MDS2, color = Month, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS1")) +
  ylab(paste("MDS2")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_MONTH_k3_MDS12

ggsave("EF2022NMDS_MONTH_k3_MDS12.tiff")
## Saving 7 x 5 in image
#MDS1 vs. MDS3:
EF2022NMDS_MONTH_k3_MDS13 = ggplot(site.scores, aes(x=MDS1, y=MDS3, color = Month, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS1")) +
  ylab(paste("MDS3")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_MONTH_k3_MDS13

ggsave("EF2022NMDS_MONTH_k3_MDS13.tiff")
## Saving 7 x 5 in image
#MDS2 vs. MDS3
EF2022NMDS_MONTH_k3_MDS23 = ggplot(site.scores, aes(x=MDS2, y=MDS3, color = Month, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS2")) +
  ylab(paste("MDS3")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_MONTH_k3_MDS23

ggsave("EF2022NMDS_MONTH_k3_MDS23.tiff")
## Saving 7 x 5 in image
## VEGETATION ZONE
# Color code by MONTH, code shapes by Unit.

# MDS1 vs. MDS2:
EF2022NMDS_ZONE_k3_MDS12 = ggplot(site.scores, aes(x=MDS1, y=MDS2, color = Vegetation.Type, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS1")) +
  ylab(paste("MDS2")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_ZONE_k3_MDS12

ggsave("EF2022NMDS_ZONE_k3_MDS12.tiff")
## Saving 7 x 5 in image
#MDS1 vs. MDS3:
EF2022NMDS_ZONE_k3_MDS13 = ggplot(site.scores, aes(x=MDS1, y=MDS3, color = Vegetation.Type, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS1")) +
  ylab(paste("MDS3")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_ZONE_k3_MDS13

ggsave("EF2022NMDS_ZONE_k3_MDS13.tiff")
## Saving 7 x 5 in image
#MDS2 vs. MDS3
EF2022NMDS_ZONE_k3_MDS23 = ggplot(site.scores, aes(x=MDS2, y=MDS3, color = Vegetation.Type, shape = Unit)) +
  geom_point(size = 3) +
  #coord_cartesian(xlim = c(-0.25,0.8), ylim = c(-0.25,0.8)) +
  scale_color_viridis_d() +
  scale_shape_manual(values = c(15,16,17,18,19)) +
  #stat_ellipse(level = 0.50) +
  theme_bw() +
  xlab(paste("MDS2")) +
  ylab(paste("MDS3")) +
  ggtitle("Electrofishing NMDS plot using Bray-Curtis Dissimilarity")
  #annotate(geom="text", x=-0.38, y=0.62, label="(k = 3, stress = 0.145)", color="black")

EF2022NMDS_ZONE_k3_MDS23

ggsave("EF2022NMDS_ZONE_k3_MDS23.tiff")
## Saving 7 x 5 in image
# ADONIS (aka. PERMANOVA) was chosen as significance test due to (1) our data are multivariate,(2) our data are nonparametric, and (3) the distance matrix we created uses Bray-Curtis distance.

EF2022ADONIS = adonis2((EF_distmat) ~ Unit + Vegetation.Type + Month + Unit:Month + Unit:Vegetation.Type, data = site.scores, permutations = 999, method = "bray")
summary(EF2022ADONIS)
##        Df            SumOfSqs             R2                F        
##  Min.   : 1.000   Min.   : 0.1459   Min.   :0.01398   Min.   :0.691  
##  1st Qu.: 2.000   1st Qu.: 0.5475   1st Qu.:0.05248   1st Qu.:1.115  
##  Median : 4.000   Median : 2.5334   Median :0.24285   Median :1.479  
##  Mean   : 8.286   Mean   : 2.9806   Mean   :0.28571   Mean   :1.930  
##  3rd Qu.:10.000   3rd Qu.: 3.3289   3rd Qu.:0.31910   3rd Qu.:1.520  
##  Max.   :29.000   Max.   :10.4320   Max.   :1.00000   Max.   :4.845  
##                                                       NA's   :2      
##      Pr(>F)      
##  Min.   :0.0010  
##  1st Qu.:0.0250  
##  Median :0.0690  
##  Mean   :0.2322  
##  3rd Qu.:0.3280  
##  Max.   :0.7380  
##  NA's   :2
print(EF2022ADONIS)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = (EF_distmat) ~ Unit + Vegetation.Type + Month + Unit:Month + Unit:Vegetation.Type, data = site.scores, permutations = 999, method = "bray")
##                      Df SumOfSqs      R2      F Pr(>F)    
## Unit                  4   4.0910 0.39216 4.8445  0.001 ***
## Vegetation.Type       1   0.1459 0.01398 0.6910  0.738    
## Month                 2   0.4707 0.04512 1.1148  0.328    
## Unit:Month            8   2.5668 0.24605 1.5198  0.025 *  
## Unit:Vegetation.Type  2   0.6243 0.05985 1.4786  0.069 .  
## Residual             12   2.5334 0.24285                  
## Total                29  10.4320 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1