Veg_Ordination_Methods_Analysis
Eli Weaver, Alexandra Bozimowski
2023/1/26
Vegetation NMDS Analysis for 2023 Shiawassee Master’s Project
Team
University of Michigan School for Environment and Sustainability
Loading packages
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.2library(vegan)## Loading required package: permute## Warning: package 'permute' was built under R version 4.1.2## Loading required package: lattice## This is vegan 2.5-7library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.2## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──## ✔ tibble 3.1.6 ✔ dplyr 1.0.10
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ✔ purrr 1.0.1## Warning: package 'tidyr' was built under R version 4.1.2## Warning: package 'readr' was built under R version 4.1.2## Warning: package 'purrr' was built under R version 4.1.2## Warning: package 'dplyr' was built under R version 4.1.2## Warning: package 'stringr' was built under R version 4.1.2## Warning: package 'forcats' was built under R version 4.1.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()Prepare data.
# Read in an original data file with IVI values
# Columns: unit (transform from 'Location'), veg (transform 'Vegetation Type'), species (transform 'Plant Name'), IVI
VegNMDSprep = read.csv("2023_official_veg_IVI.csv", header = T)
# Create a 'siteID' column for unique observations and remove NA values
VegNMDSprep.2 = VegNMDSprep%>%
unite(siteID, unit, veg, remove = F)%>%
filter(!is.na(IVI))
# Create your working NMDS data frame which uses the 'siteID' column you created above
IVInmds = pivot_wider(
VegNMDSprep.2,
id_cols = siteID,
names_from = species,
names_repair = "check_unique",
values_from = IVI,
values_fill = 0)
IVInmds## # A tibble: 17 × 93
## siteID `Abutilon theo…` `Alisma planta…` `Amaranthus sp…` `Ammannia cocc…`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 MC_ME 5.30 2.16 1.07 1.07
## 2 MC_Phala… 0 0 0 0
## 3 MC_Salix 0 0 0 0
## 4 MC_SAV 0 0 0 0
## 5 MC_Typha 0 0 0 0
## 6 MN_Phala… 14.8 2.97 0 0
## 7 MN_SAV 0 0 0 0
## 8 MN_Typha 11.8 2.41 1.98 0
## 9 MS_Fores… 0 0 0 0
## 10 MS_Fores… 0 6.30 0 0
## 11 MS_ME 0 3.19 0 0
## 12 MS_SAV 0 0 0 0
## 13 MS_Typha 0 0 0 0
## 14 P1A_Mudf… 19.1 16.4 0 0
## 15 P1A_Nymp… 0 4.23 0 0
## 16 P1A_Salix 1.72 0 0 0
## 17 P1A_Typha 0 0 0 0
## # … with 88 more variables: `Bidens cernua` <dbl>,
## # `Ceratophyllum demersum` <dbl>, `Cirsium arvense` <dbl>,
## # `Cyperus odorata` <dbl>, `Dead Typha` <dbl>, `Eleocharis acicularis` <dbl>,
## # `Eleocharis palustris` <dbl>, `Glebionis coronaria` <dbl>,
## # `Lemna minor` <dbl>, `Lemna trisulca` <dbl>, Moss <dbl>,
## # `Najas minor` <dbl>, `Oenothera laciniata` <dbl>, `Panicum spp.` <dbl>,
## # `Persicaria lapathifolia` <dbl>, `Persicaria maculosa` <dbl>, …# Write this file back to your working directory so that you can create your grouping spreadsheet(s)
write.csv(IVInmds, "UMSEAS-SNWR-2022_Vegetation-NMDS-PERMANOVA.csv")
# Reopen "UMSEAS-SNWR-2022_Vegetation-NMDS-PERMANOVA.csv," delete the IVI values columns and ADD a 'unit' column and a 'veg' column.
# Save as "UMSEAS-SNWR-2022_veg_NMDS_trts.csv";
veg_trts = read.csv("UMSEAS-SNWR-2022_veg_NMDS_trts.csv", header = T)
# Remove the 'siteID' column to perform the next steps
IVInmds = IVInmds[,2:93]# Calculate user-defined distance matrix using "method=Bray-Curtis Dissimilarity".
# Note: Result matrix from this block will be later used in NMDS.
distmat <- vegdist(IVInmds, method = "bray")
distance.matrix <- as.matrix(distmat, labels=T)
#distance.matrix ## Result is too long for the knitted html. Plz uncomment this line of code if you want.
# Write this file back to your working directory
write.csv(distance.matrix, "UMSEAS-SNWR-2022_Veg_bray_distmatrix.csv")5) Draw NMDS plot(s) using the full IVI Matrix and the Bray-Curtis Dissimilarity Matrix.
# 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
vegNMDS <- metaMDS(distance.matrix, distance="bray",k=3,maxit=999,trymax=500,wascores=TRUE)## Run 0 stress 0.07293249
## Run 1 stress 0.07546485
## Run 2 stress 0.07871289
## Run 3 stress 0.07546562
## Run 4 stress 0.07293307
## ... Procrustes: rmse 0.0003245438 max resid 0.0008427241
## ... Similar to previous best
## Run 5 stress 0.07293264
## ... Procrustes: rmse 0.001906031 max resid 0.004816333
## ... Similar to previous best
## Run 6 stress 0.07871091
## Run 7 stress 0.08252295
## Run 8 stress 0.07871259
## Run 9 stress 0.07871227
## Run 10 stress 0.08252352
## Run 11 stress 0.07871291
## Run 12 stress 0.07293332
## ... Procrustes: rmse 0.00288687 max resid 0.007784385
## ... Similar to previous best
## Run 13 stress 0.07293367
## ... Procrustes: rmse 0.002556258 max resid 0.007543925
## ... Similar to previous best
## Run 14 stress 0.07546465
## Run 15 stress 0.07293262
## ... Procrustes: rmse 0.002422384 max resid 0.006513866
## ... Similar to previous best
## Run 16 stress 0.0754656
## Run 17 stress 0.07871318
## Run 18 stress 0.07871295
## Run 19 stress 0.07293297
## ... Procrustes: rmse 0.001956506 max resid 0.005834165
## ... Similar to previous best
## Run 20 stress 0.07871219
## *** Solution reachedgoodness(vegNMDS)## [1] 0.009797266 0.016188403 0.021124228 0.012603738 0.026527469 0.015782629
## [7] 0.008787387 0.012026372 0.010213539 0.025397407 0.019655876 0.018937533
## [13] 0.008129121 0.023632637 0.016960704 0.023993783 0.014109587stressplot(vegNMDS)
vegNMDS##
## Call:
## metaMDS(comm = distance.matrix, distance = "bray", k = 3, trymax = 500, wascores = TRUE, maxit = 999)
##
## global Multidimensional Scaling using monoMDS
##
## Data: distance.matrix
## Distance: user supplied
##
## Dimensions: 3
## Stress: 0.07293249
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation
## Species: scores missing# Create a new df that gets your MDS output data for plotting by adding the 'points' column to the vegNMDS df.
site.scores = as_tibble(vegNMDS$points)
site.scores## # A tibble: 17 × 3
## MDS1 MDS2 MDS3
## <dbl> <dbl> <dbl>
## 1 -0.0215 0.0658 0.0229
## 2 0.224 0.128 0.0566
## 3 0.0408 -0.102 -0.0866
## 4 -0.197 -0.210 0.0281
## 5 -0.139 0.108 -0.107
## 6 0.134 0.236 0.156
## 7 -0.234 -0.208 0.0331
## 8 -0.173 0.290 0.0560
## 9 0.705 -0.186 -0.00457
## 10 0.143 -0.0149 -0.00801
## 11 -0.143 -0.0971 -0.122
## 12 -0.152 -0.166 -0.103
## 13 -0.224 0.0395 -0.0447
## 14 -0.0847 -0.0210 0.271
## 15 -0.0403 -0.219 0.115
## 16 0.247 0.129 -0.166
## 17 -0.0858 0.229 -0.0976# 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(veg_trts$unit))
site.scores$Zone = factor(paste(veg_trts$veg))
# Time to plot!
# Color code by Vegetation Zone, code shapes by Unit.
# MDS1 vs. MDS2: (The first time, run without 'coord_cartesian()' and 'annotate()', then add those customizations 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.)
Veg2022NMDS_k3_MDS12 = ggplot(site.scores, aes(x=MDS1, y=MDS2, color=Zone, shape=Unit))+
geom_point(size = 3)+
coord_cartesian(xlim = c(-0.6,0.4), ylim = c(-0.6,0.4))+
scale_color_viridis_d()+
scale_shape_manual(values = c(15,16,17,18))+
stat_ellipse(level = 0.50)+
theme_bw()+
xlab(paste("MDS1")) +
ylab(paste("MDS2")) +
ggtitle("Vegetation NMDS plot using Bray-Curtis Dissimilarity")
#geom_text_repel(aes(label = paste0(Unit, " ", Zone), fontface = 2), nudge_y = 0.04, size = 3) +
#annotate(geom="text", x=-0.48, y=0.42, label="(k = 3, stress = 0.094)", color="black")
Veg2022NMDS_k3_MDS12
ggsave("UMSEAS-SNWR-2022_vegNMDS_k3_MDS12.tiff")
#MDS1 vs. MDS3:
Veg2022NMDS_k3_MDS13 = ggplot(site.scores, aes(x=MDS1, y=MDS3, color=Zone, shape=Unit))+
geom_point(size = 3)+
coord_cartesian(xlim = c(-0.6,0.4), ylim = c(-0.6,0.4))+
scale_color_viridis_d()+
scale_shape_manual(values = c(15,16,17,18))+
stat_ellipse(level = 0.50)+
theme_bw()+
xlab(paste("MDS1")) +
ylab(paste("MDS3")) +
ggtitle("Vegetation NMDS plot using Bray-Curtis Dissimilarity")
#geom_text_repel(aes(label = paste0(Unit, " ", Zone), fontface = 2), nudge_y = 0.04, size = 3) +
#annotate(geom="text", x=-0.48, y=0.42, label="(k = 3, stress = 0.094)", color="black")
Veg2022NMDS_k3_MDS13
ggsave("UMSEAS-SNWR-2022_vegNMDS_k3_MDS13.tiff")
#MDS2 vs. MDS3:
Veg2022NMDS_k3_MDS23 = ggplot(site.scores, aes(x=MDS2, y=MDS3, color=Zone, shape=Unit))+
geom_point(size = 3)+
coord_cartesian(xlim = c(-0.6,0.4), ylim = c(-0.6,0.4))+
scale_color_viridis_d()+
scale_shape_manual(values = c(15,16,17,18))+
stat_ellipse(level = 0.50)+
theme_bw()+
xlab(paste("MDS2")) +
ylab(paste("MDS3")) +
ggtitle("Vegetation NMDS plot using Bray-Curtis Dissimilarity")
#geom_text_repel(aes(label = paste0(Unit, " ", Zone), fontface = 2), nudge_y = 0.04, size = 3) +
#annotate(geom="text", x=-0.48, y=0.42, label="(k = 3, stress = 0.092)", color="black")
Veg2022NMDS_k3_MDS23
ggsave("UMSEAS-SNWR-2022_vegNMDS_k3_MDS23.tiff")
#Graph 1v2 is most representative of significant difference in IVI values based upon vegetation zones. Color coded veg zones clearly grouped together, whereas we did not see any significance by unit and thus do not see grouping based on shape. ##5.1) Dealing with (NMDS) STRESS. In the event asking the NMDS to resolve to 3 dimensions (i.e., k=3) results in excessive stress, you can modify the resolution. For example, in 2021 they tested k=2. Three is best. Four gets excessive but is doable. Two is a bit silly.
# Below are code for the case "k=2", because: (1) when k=2, stress = 0.151 which is fair; (2) easy to interpret in 2D plot
#vegNMDS_2 <- metaMDS(distance.matrix, distance="bray",k=2,maxit=999,trymax=500,wascores=TRUE)
#goodness(vegNMDS_2)
#stressplot(vegNMDS_2)
#site.scores_2 = as_tibble(vegNMDS_2$points)
#site.scores_2$Unit = factor(paste(veg_trts$unit), labels = c("MC", "MN", "MS", "P1A"))
#site.scores_2$Zone = factor(paste(veg_trts$zone), labels = c("MEV", "Phalaris", "Salix", "SAV", "Typha", "Forest", "Nymphaea"))
#MDS1 vs. MDS2:
#Veg2021NMDS_k2 = ggplot(site.scores_2, aes(x=MDS1, y=MDS2, color=Zone, shape=Unit)) +
#geom_point(size = 3) +
#coord_cartesian(xlim = c(-0.6,0.4), ylim = c(-0.6,0.4)) +
#scale_color_viridis_d() +
#scale_shape_manual(values = c(15,16,17,18)) +
#stat_ellipse(level = 0.50) +
#theme_bw() +
#xlab(paste("MDS1")) +
#ylab(paste("MDS2")) +
#ggtitle("Vegetation NMDS plot using Bray-Curtis Dissimilarity") +
#geom_text_repel(aes(label = paste0(Unit, " ", Zone), fontface = 2), nudge_y = 0.04, size = 3) +
#annotate(geom="text", x=-0.48, y=0.42, label="(k = 2, stress = 0.151)", color="black")
#Veg2021NMDS_k2
#ggsave("UMSEAS-SNWR-2021_vegNMDS_k2.tiff")# 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.
#Veg2022ADONIS = adonis2((distance.matrix) ~ Unit + Zone + Unit:Zone, data = site.scores, permutations = 999, method = "bray")
#summary(Veg2022ADONIS)
#print(Veg2022ADONIS)
#Sample size was too low to run an interaction test, ie unit;zone portion of the adonis testing. More samples would be needed in order to run this effectively
Veg2022ADONIS_unit = adonis2((distance.matrix) ~ Unit, data = site.scores, permutations = 999, method = "bray")
summary(Veg2022ADONIS_unit)## Df SumOfSqs R2 F
## Min. : 3.00 Min. :0.9473 Min. :0.1808 Min. :0.9564
## 1st Qu.: 8.00 1st Qu.:2.6197 1st Qu.:0.5000 1st Qu.:0.9564
## Median :13.00 Median :4.2921 Median :0.8192 Median :0.9564
## Mean :10.67 Mean :3.4930 Mean :0.6667 Mean :0.9564
## 3rd Qu.:14.50 3rd Qu.:4.7658 3rd Qu.:0.9096 3rd Qu.:0.9564
## Max. :16.00 Max. :5.2394 Max. :1.0000 Max. :0.9564
## NA's :2
## Pr(>F)
## Min. :0.546
## 1st Qu.:0.546
## Median :0.546
## Mean :0.546
## 3rd Qu.:0.546
## Max. :0.546
## NA's :2print(Veg2022ADONIS_unit)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = (distance.matrix) ~ Unit, data = site.scores, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Unit 3 0.9473 0.1808 0.9564 0.546
## Residual 13 4.2921 0.8192
## Total 16 5.2394 1.0000Veg2022ADONIS_zone = adonis2((distance.matrix) ~ Zone, data = site.scores, permutations = 999, method = "bray")
summary(Veg2022ADONIS_zone)## Df SumOfSqs R2 F
## Min. : 8.00 Min. :1.209 Min. :0.2307 Min. :3.335
## 1st Qu.: 8.00 1st Qu.:2.620 1st Qu.:0.5000 1st Qu.:3.335
## Median : 8.00 Median :4.031 Median :0.7693 Median :3.335
## Mean :10.67 Mean :3.493 Mean :0.6667 Mean :3.335
## 3rd Qu.:12.00 3rd Qu.:4.635 3rd Qu.:0.8847 3rd Qu.:3.335
## Max. :16.00 Max. :5.239 Max. :1.0000 Max. :3.335
## NA's :2
## Pr(>F)
## Min. :0.001
## 1st Qu.:0.001
## Median :0.001
## Mean :0.001
## 3rd Qu.:0.001
## Max. :0.001
## NA's :2print(Veg2022ADONIS_zone)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = (distance.matrix) ~ Zone, data = site.scores, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Zone 8 4.0308 0.76932 3.3351 0.001 ***
## Residual 8 1.2086 0.23068
## Total 16 5.2394 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1