library(tidyverse)   
library(readxl)
library(rstatix)
library(fitdistrplus) 
library(lme4) 
library(bbmle) 
library(RVAideMemoire) 
library(car) 
library(emmeans)
library(multcomp) 
library(ggpubr)
library(plotly)
library(gridExtra)
library(grid)


#Latency to target area - sheet 2
Latency.Maze <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 2)

attach(Latency.Maze)

descdist(Latency, discrete=FALSE) 
plot(density(Latency))
hist(Latency)

######
#Fig 3a

LaTr1 <- glmer(Latency ~ Group * Day + (1 | ID),
            data = Latency.Maze,
            family=Gamma(link = "log"))

LaTr2 <- glmer(Latency ~ Group + Day + (1 | ID),
            data = Latency.Maze,
            family=Gamma(link = "log"))

AICctab(LaTr1, LaTr2)
plotresid(LaTr2, shapiro = T)
summary(LaTr2)
Anova(LaTr2, test.statistic="Chisq", type="III")
#poshoc 
ph.LT <- emmeans(LaTr2, ~ Group, adjust="bonferroni")
pairs(ph.LT)
ph.LT2 <- emmeans(LaTr2, ~ Day, adjust="bonferroni")
pairs(ph.LT2)

Latency.Maze$Group <- factor(Latency.Maze$Group, levels = c("Control", "Moderate warming", "High warming"))

#Parallel plot or parallel coordinates plot allows to compare the feature of several individual observations (series) on a set of numeric variables. Each vertical bar represents a variable and often has its own scale. (The units can even be different). Values are then plotted as series of lines connected across each axis.
library(GGally)
Fig3a <- ggparcoord(Latency.Maze, scale = "uniminmax",
                    columns = c(2, 5),
                    groupColumn = 1,
                    order = "allClass",
                    alphaLines = 0.1, boxplot = TRUE,
                    showPoints = F,
                    title="(a) associative learning") +
  facet_wrap(~ Group) + coord_flip() +
  scale_x_discrete(name ="Fixed effects") +
  scale_color_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) +
  theme_bw()
Fig3a


######
#Fig 3b
datLearMem <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 3)

descdist(datLearMem$Latency, discrete=FALSE) 
plot(density(datLearMem$Latency))
hist(datLearMem$Latency)

LearMem1 <- glmer(Latency ~ Group * Day + (1 | ID),
               data = datLearMem,
               family=Gamma(link = "log"))

LearMem2 <- glmer(Latency ~ Group + Day + (1 | ID),
               data = datLearMem,
               family=Gamma(link = "log"))

AICctab(LearMem1, LearMem2)
plotresid(LearMem2, shapiro = T)
summary(LearMem2)
Anova(LearMem2, test.statistic="Chisq", type="III")
#poshoc 
ph.LM <- emmeans(LearMem2, ~ Group, adjust="bonferroni")
pairs(ph.LM)
ph.LM2 <- emmeans(LearMem2, ~ Day, adjust="bonferroni")
pairs(ph.LM2)

datLearMem$Group <- factor(datLearMem$Group, levels = c("Control", "Moderate warming", "High warming"))

Fig3b <- ggpaired(datLearMem, x = "Day", y = "Latency",
         fill = "Group", color = "black",
         legend = "none", xlab ="",
         title = "(b) learning and memory test",
         ylab = "Latency to enter the target hall (s)",
         palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
         line.color = "#5b5b5b", line.size = 0.4,
         short.panel.labs = T) +
         facet_wrap(~Group, ncol = 3)
Fig3b 

ggarrange(Fig3a, Fig3b, ncol = 1, nrow = 2,
          common.legend = T, legend = "bottom")