# Load libraries
library(dplyr)
library(tidyr)
library(broom)
library(ggpubr)
library(emmeans)
library(ordinal)
# install.packages("viridis")  # Install para los colors
#library("viridis")           # Load



# TORTULA as Target (The same must be run for Syntrichia file)

# Load Data
nc <- read.csv("file:///C:/Users/Data/T_target.csv", sep=";")

nc$plate <- rownames(nc)

molten <- gather(nc, Class, value, -Culture, -Emmiter, -Target, -Treatment,-plate)
molten$value <- as.numeric(molten$value)
molten_nc <- molten[rep(seq_len(nrow(molten)), 
                           molten$value), ][, c(1:6)]

molten_nc$sizeMedian = with(molten_nc, 
                         ifelse(Class == "Unicelular", 1,
                                ifelse(Class == "P2_20", 11,
                                       ifelse(Class == "P21_40", 30,
                                              ifelse(Class == "P41_100", 70, 
                                                     ifelse(Class == "P100", 300, "NA")))))
)

molten_nc$sizeMedian <- as.numeric(molten_nc$sizeMedian)

# Filter the data, from here we need to repeat the code for every experiment
es_to   <- dplyr::filter(molten_nc, Target=="To"& Emmiter=="Es")
# Take out the leves that are NaN
es_to   <- droplevels(es_to  )
# Order the levels
es_to  $Treatment <- factor(es_to $Treatment, levels=c("Control","Es_Sy","Es_To"))
# Set the categories order
es_to  $Class <- factor(es_to  $Class, ordered = TRUE, levels = c("Unicelular","P2_20","P21_40","P41_100", "P100"))
# Transform into a factor the variable "plate" (random effect for every Nalgen Jar)
es_to  $plate <- as.factor(es_to  $plate)

# Calculate the summary of spores per treatment
ftable(xtabs(~ Target + Emmiter + Treatment + Class, data = es_to  ))

# ANOVA
m1 <- nlme::lme(sizeMedian ~ Treatment, random=~ 1|plate, data = es_to  )
m1
anova(m1)

# posthoc
posthoc <- emmeans(m1, pairwise~Treatment, adjust="tukey")
tposthoc <- tidy(posthoc$contrasts)
posthoc
write.csv(tposthoc, "es_to  _posth.csv")
# getwd() 

# ordered logistic regression

fm1 <- clmm(Class ~ Treatment + (1|es_to  $plate), data = es_to  )
fm1
summary(fm1)
emmeans <- emmeans(fm1, "Treatment", type = "response")
emmeans
pairs <- pairs(emmeans, reverse=TRUE)
pairs
write.csv(pairs, "logposh.csv")