###Data frame: raw abundance of periphyton.
library(readr)
library(dplyr)
library(tidyverse)
library(tibble)
praw <- read.csv("PuntosPeriRaw.csv")
praw$Land.use <- as.factor(praw$Land.use)
praw$River <- as.factor(praw$River)
praw$Point <- as.factor(praw$Point)
for (i in 4:ncol(praw)) {praw[[i]] <- as.numeric(praw[[i]])}
str(praw)
##Subset
peri <- praw %>%
  dplyr::select(-Land.use, -River) %>%
  column_to_rownames("Point")
peri
site <- praw %>%
  dplyr::select(Land.use, River, Point)
site

###Biodiversity information
library(iNEXT)
peri <- t(peri)
##Calculating S, H, D and J with iNEXT
S <- ChaoRichness(peri, datatype = "abundance", conf = 0.95)
S <- S %>%
  rename(S.Observed = Observed,
         S.Estimator = Estimator,
         S.Est_s.e = Est_s.e.,
         S.Lower = `95% Lower`,
         S.Upper = `95% Upper`)

H <- ChaoShannon(peri, datatype = "abundance", transform = FALSE, conf = 0.95, B = 200)
H <- H %>%
  rename(H.Observed = Observed,
         H.Estimator = Estimator,
         H.Est_s.e = Est_s.e,
         H.Lower = `95% Lower`,
         H.Upper = `95% Upper`)

D <- ChaoSimpson(peri, datatype = "abundance", transform = FALSE, conf = 0.95, B = 200)
D <- D %>%
  rename(D.Observed = Observed,
         D.Estimator = Estimator,
         D.Est_s.e = Est_s.e.,
         D.Lower = `95% Lower`,
         D.Upper = `95% Upper`)

J <- (H$H.Estimator/log(S$S.Estimator))

bio <- cbind(S, H, D, J)
bio <- as.data.frame(bio)
bio <- rownames_to_column(bio, var = "Point")
bio <- dplyr::left_join(site, bio, by = "Point")
bio
write.csv(bio, "iNEXT_biodiversityindex.csv")

##Comparisons between land uses
#lme, anova and Tukey test
library(nlme)
library(multcomp)
#Richness
S.a <- lme(S.Estimator ~ Land.use, data = bio, random = ~ 1 | River)
summary(S.a)
qqnorm(resid(S.a))
qqline(resid(S.a))
anova(S.a) #p-value land use 0.0078#
local({
  .Pairs <- glht(S.a, linfct = mcp(Land.use = "Tukey"))
  print(summary(.Pairs)) # pairwise tests
  print(confint(.Pairs, level=0.95)) 
  print(cld(.Pairs, level=0.05)) 
  old.oma <- par(oma=c(0, 5, 0, 0))
  plot(confint(.Pairs))
  par(old.oma)
}) #BOs(a), GLs(b), OPs(ab), OPsZr(a)#
ggplot(bio, aes(x=Land.use, y=S.Estimator, fill= Land.use)) + 
  geom_boxplot() + 
  labs(x = NULL, y = "Richness", fill = "Land use") +
  scale_fill_manual(values=c("#088908", "#980d0d", "#e98e90", "yellow3")) +
  theme_minimal() + theme(legend.position = "bottom")

#Simpson
D.a <- lme(D.Estimator ~ Land.use, data = bio, random = ~ 1 | River)
summary(D.a)
qqnorm(resid(D.a))
qqline(resid(D.a))
anova(D.a) #p-value land use 0.0423
local({
  .Pairs <- glht(D.a, linfct = mcp(Land.use = "Tukey"))
  print(summary(.Pairs)) # pairwise tests
  print(confint(.Pairs, level=0.95)) 
  print(cld(.Pairs, level=0.05)) 
  old.oma <- par(oma=c(0, 5, 0, 0))
  plot(confint(.Pairs))
  par(old.oma)
}) #BOs(ab), GLs(b), OPs(a), OPsZr(ab)#
ggplot(bio, aes(x=Land.use, y=D.Observed, fill= Land.use)) + 
  geom_boxplot() + 
  labs(x = NULL, y = "Simpson Index", fill = "Land use") +
  scale_fill_manual(values=c("#088908", "#980d0d", "#e98e90", "yellow3")) +
  theme_minimal() + theme(legend.position = "bottom")

#Pielou
J.a <- lme(J ~ Land.use, data = bio, random = ~ 1 | River)
summary(J.a)
qqnorm(resid(J.a))
qqline(resid(J.a))
anova(J.a) #p-value land use 0.1883
ggplot(bio, aes(x=Land.use, y=J, fill= Land.use)) + 
  geom_boxplot() + 
  labs(x = NULL, y = "Pielou Evenness Index", fill = "Land use") +
  scale_fill_manual(values=c("#088908", "yellow3", "#980d0d", "#e98e90")) +
  theme_minimal() + theme(legend.position = "bottom")
min(bio$J) #0.4441318
max(bio$J) #0.9423188

#Visualization
barp <- read.csv("barplot_biod.csv")
barp$Land.use <- as.factor(barp$Land.use)
barp$Index <- as.factor(barp$Index)
bar <- barp %>%
  slice(1:4,)
bar
ggplot(barp, aes(x = Land.use, y = Mean, fill = Land.use)) + 
  geom_bar(stat = "identity", position = position_dodge()) +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), 
                width = 0.2, position = position_dodge(0.9)) +
  scale_fill_manual(values = c("FO" = "#088908", "PA" = "yellow3", "OP" = "#980d0d", "OPB" = "#e98e90")) +
  guides(fill = guide_legend(override.aes = list(color = NULL))) +
  labs(color = NULL) + 
  facet_wrap(~Index, scales = "free", ncol = 3, labeller = label_value) +
  theme_minimal() + theme(legend.position = "bottom") + 
  labs(x = NULL, y = NULL, fill = "Land use")