nMDS TRCA 2011

library(devtools)
library(rmarkdown)
library(vegan)
library(tidyverse)
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(ggalt)
library(knitr)

Import

TRCA_base <- read.csv("TRCA_Herbaceous.csv")
CLOCA_base <- read.csv("CLOCA_base.csv")
CVC_base <- read.csv("CVC_base.csv")

Cleaning

# Cleaning of TRCA
TRCA_clean <- TRCA_base %>% 
  select(SiteName, PlotID, Subplot, Date, LandscapeCategory, Species, CommonName, PlantType, Cover) %>% 
  mutate(Cover = as.numeric(ifelse(Cover == "<1", 0.5, Cover))) %>% 
  mutate(Plot_ent = paste(PlotID, Subplot, sep = "-sub"))


# Cleaned dataset, subset for 2011
TRCA_2011 <- TRCA_clean %>% 
  filter(year(as.Date(Date, format = "%d-%b-%y")) == 2011) %>% 
  select(Species, Cover, Plot_ent, SiteName, PlotID, LandscapeCategory) %>%
  group_by(Species, Plot_ent) %>% 
  mutate(mean_value = ifelse(length(Cover) > 1, mean(Cover), Cover)) %>% 
  distinct(Species, Plot_ent, .keep_all = TRUE) %>% 
  select(-Cover) %>% 
  rename(Cover = mean_value)

2011 nMDS

# Pivoting 2011 data to be wide
# Creates unique identifier for subplot (Plot_ent) var
# Replaces NA with 0
pivot_TRCA_2011 <- TRCA_2011 %>% 
  pivot_wider(names_from = Species, values_from = Cover) %>% 
  filter(!duplicated(Plot_ent)) %>% 
  column_to_rownames(var = "Plot_ent") %>%
  replace(is.na(.), 0.0)


# Selects for the species data within pivot
matrix_TRCA_2011 <- pivot_TRCA_2011 %>% 
  select(4:284)


# nMDS function
NMDS_TRCA_2011 <- metaMDS(matrix_TRCA_2011, 
                          distance = "bray", k = 2, autotransform = FALSE)

## Run 0 stress 0.1324329 
## Run 1 stress 0.1361502 
## Run 2 stress 0.1334549 
## Run 3 stress 0.1352116 
## Run 4 stress 0.1354962 
## Run 5 stress 0.133071 
## Run 6 stress 0.1358424 
## Run 7 stress 0.1340179 
## Run 8 stress 0.1353289 
## Run 9 stress 0.1340589 
## Run 10 stress 0.1365971 
## Run 11 stress 0.1347093 
## Run 12 stress 0.1346124 
## Run 13 stress 0.134436 
## Run 14 stress 0.1359184 
## Run 15 stress 0.132882 
## ... Procrustes: rmse 0.01853122  max resid 0.103534 
## Run 16 stress 0.14038 
## Run 17 stress 0.1352785 
## Run 18 stress 0.1321213 
## ... New best solution
## ... Procrustes: rmse 0.02247795  max resid 0.1248595 
## Run 19 stress 0.1351005 
## Run 20 stress 0.1345196 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax

# Combining nMDS points with ecological data
plot_df_NMDS_TRCA_2011 <- data.frame(
  NMDS_TRCA_2011$points,
  Landscape = pivot_TRCA_2011$LandscapeCategory,
  Site = pivot_TRCA_2011$SiteName)


# Plotting of nMDS
ggplot(plot_df_NMDS_TRCA_2011, 
       aes(x = MDS1, y = MDS2, color = Landscape)) +
  geom_point(size = 1) +
  coord_fixed() +
  scale_color_manual(values = c("blue", "red", "green")) +
  theme_classic() +
  theme(legend.position = "right") +
# Circles points based on site
  stat_ellipse(aes(x = MDS1, y = MDS2, color = Landscape), 
                alpha = 0.2) +
# Prints stress onto plot
  annotate(geom = "label", x = -2, y = 1, size = 5,
           label = paste("Stress: ", 
                         round(NMDS_TRCA_2011$stress, digits = 3))) +
# Print title
  labs(title = "nMDS TRCA 2011")

# Creates Shepard plot
stressplot(NMDS_TRCA_2011, main = "Stress plot TRCA 2011")

# Analysis of variance using distance matrices
fit_NMDS_TRCA_2011 <- adonis2(matrix_TRCA_2011 ~ LandscapeCategory, data = pivot_TRCA_2011, permutations = 999, method = "bray")

fit_NMDS_TRCA_2011 

# Bray-Curtis distances between samples
distances_TRCA_2011 <- vegdist(matrix_TRCA_2011)

disper_TRCA_2011 <- betadisper(distances_TRCA_2011, pivot_TRCA_2011$LandscapeCategory)

disper_TRCA_2011

## 
##  Homogeneity of multivariate dispersions
## 
## Call: betadisper(d = distances_TRCA_2011, group =
## pivot_TRCA_2011$LandscapeCategory)
## 
## No. of Positive Eigenvalues: 145
## No. of Negative Eigenvalues: 93
## 
## Average distance to median:
##      rural      urban urbanizing 
##     0.6567     0.6444     0.5810 
## 
## Eigenvalues for PCoA axes:
## (Showing 8 of 238 eigenvalues)
##  PCoA1  PCoA2  PCoA3  PCoA4  PCoA5  PCoA6  PCoA7  PCoA8 
## 16.649  9.973  7.269  5.756  3.990  3.756  3.492  3.278

anova(disper_TRCA_2011)