###Data frame: chlorophyll mg/m2
library(readr)
library(dplyr)
library(tidyverse)
#Visualization
library(ggplot2)
chlap <- read.csv("Chlorophyll.plot.csv", na.strings = "NA")
chlap$Land.use <- as.factor(chlap$Land.use)
chlap$River <- as.factor(chlap$River)
chlap$Sample.type <- as.factor(chlap$Sample.type)
chlap$Chlorophyl <- as.numeric(chlap$Chlorophyl)
str(chlap)
ggplot(chlap, aes(x=Land.use, y=Chlorophyl, fill= Land.use)) + 
  geom_boxplot() + 
  geom_jitter(shape=16, position=position_jitter(0.2), alpha = 0.5) +
  scale_fill_manual(values=c("FO" = "#088908", "PA" = "yellow3", "OP" = "#980d0d", "OPB" = "#e98e90")) + 
  facet_wrap(~Sample.type, scales = "free", ncol = 2) +
  theme_minimal() + theme(legend.position = "bottom") + 
  labs(x = NULL, y = "Chlorophyll mg/m2", fill = "Land use")

#Stadistical difference
cloro <- read.csv("Chlorophyll.dif.csv", na.strings = "NA")
cloro$Land.use <- as.factor(cloro$Land.use)
cloro$River <- as.factor(cloro$River)
str(cloro)
library(nlme)
library(multcomp)

Sediments <- lme(Sediment ~ Land.use, data = cloro, random = ~ 1 | River)
summary(Sediments)
qqnorm(resid(Sediments))
qqline(resid(Sediments))
anova(Sediments) #p-value land use 0.0085
local({
  .Pairs <- glht(Sediments, 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)
}) #FO(a), PA(c), OP(ab), OPB(ab)#

Stone <- lme(Stone ~ Land.use, data = cloro, random = ~ 1 | River)
summary(Stone)
qqnorm(resid(Stone))
qqline(resid(Stone))
anova(Stone) #p-value land use 0.059
local({
  .Pairs <- glht(Stone, 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)
}) #FO(ab), PA(b), OP(ab), OPB(a)#