###Data frame: density of periphyton.
library(readr)
library(dplyr)
library(tidyverse)
pabund <- read.csv("RiosPeriAbund.csv")
pabund$Land.use <- as.factor(pabund$Land.use)
pabund$Land.use1 <- as.factor(pabund$Land.use1) 
pabund$River <- as.factor(pabund$River)
for (i in 5:ncol(pabund)) {pabund[[i]] <- as.numeric(pabund[[i]])}
str(pabund)
##Subset
peri <- pabund %>%
  dplyr::select(-Land.use, -Land.use1, -Date_Colect, -DO,
         -Temperature, -Conductivity, -pH, -Turbidity,
         -DBO5, -Solids, - Silica, -PO4, -NO3, -NH4,
         -Inorganic.N, -Canopy, -Light, -ChlaSto, -ChlaSedi) %>%
  column_to_rownames("River")
site <- pabund %>% 
  dplyr::select(Land.use, Land.use1, River)
envi <- pabund %>%
  dplyr::select(River, DO, Temperature, Conductivity, pH, Turbidity,
         DBO5, Solids, Silica, NO3, NH4,
         Inorganic.N, Canopy, Light) %>%
  column_to_rownames("River")
#NMDS
library(vegan)
set.seed(123)
nmds_Peri1 <- metaMDS(peri, distance = "bray", k = 2, try = 1000, trymax = 1000)
nmds_Peri <- metaMDS(peri, distance = "bray", k = 2, try = 5000, trymax = 5000, previous.best = nmds_Peri1)
nmds_Peri$stress #0.166921
stressplot(nmds_Peri)
envfit <- envfit(nmds_Peri, envi, na = TRUE, permutations = 999) 
envfit
sppfit <- envfit(nmds_Peri, peri, permutations = 999)
sppfit
env_coord <- as.data.frame(scores(envfit, display = "vectors")) 
env_coord <- env_coord %>%
  cbind(pval = envfit$vectors$pvals) %>%
  subset(pval<=0.05) %>%
  dplyr::select(-pval)
env_coord
spp_coord <- as.data.frame(scores(sppfit, display = "vectors"))
spp_coord <- spp_coord %>%
  cbind(pval = sppfit$vectors$pvals) %>%
  subset(pval<=0.05) %>%
  dplyr::select(-pval)
spp_coord
NMDS_Comu <- as.data.frame(scores(nmds_Peri, display = "sites"))
NMDS_Comu <- NMDS_Comu %>%
  cbind("Land use" = site$Land.use, "Land use1" = site$Land.use1, River = site$River)
perinmds <- ggplot(data = NMDS_Comu, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(data = NMDS_Comu, aes(colour = `Land use`), size = 8, alpha = 1) + 
  scale_colour_manual(values = c("FO" = "#088908", "PA" = "yellow3", "OP" = "#980d0d","OPB" = "#e98e90"))  + 
  geom_text(data = env_coord, aes(x = NMDS1, y = NMDS2), colour = "grey20", 
            size= 4, fontface = "bold", label = row.names(env_coord)) +
  geom_text(aes(label = River), colour = "black", 
            size= 2, fontface = "bold") + 
  geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
               data = env_coord, size =1, alpha = 0.5, colour = "grey30") + 
  theme(axis.title = element_text(size = 10, face = "bold", colour = "grey20"), 
        panel.background = element_blank(), panel.border = element_rect(fill = NA, colour = "grey30"), 
        axis.ticks = element_blank(), axis.text = element_blank(), legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "grey30"), 
        legend.text = element_text(size = 9, colour = "grey30")) + 
  labs(colour = "Land cover")
perinmds

#PERMANOVA
perimanova <- adonis2(peri~site$Land.use, method = "bray", perm =999)
perimanova
###Pairwise permanova
install.packages("remotes")
remotes::install_github("phytomosaic/ecole")
library(ecole)
set.seed(910)
peripair1 <- permanova_pairwise(x = vegdist(peri, 'hellinger'), grp = site$Land.use1, permutations = 999, padj = "bonferroni")
peripair1
#Indicator species analysis
library(indicspecies)
set.seed(123)
spindic <-  multipatt(peri, site$Land.use1, func = "r.g", control = how(nperm=999))
summary(spindic)