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)


######################################## Fig 4a
#Fig 4a

First.Hall.Training <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 4)

attach(First.Hall.Training)


Hall1.Training <- lmer(HallNumber ~ Group * Day + (1 | ID),
               data = First.Hall.Training)

Hall2.Training <- lmer(HallNumber ~ Group + Day + (1 | ID),
               data = First.Hall.Training)

Hall3.Training <- lmer(Latency ~ HallNumber + Group + Day + (1 | ID),
              data = First.Hall.Training)

AICctab(Hall1.Training, Hall2.Training, Hall3.Training)
plotresid(Hall2.Training, shapiro = T)
summary(Hall2.Training)
Anova(Hall2.Training, test.statistic="Chisq", type="III")
#poshoc 
PostHoc.Training <- emmeans(Hall2.Training, ~ Group,
                            adjust="bonferroni")
pairs(PostHoc.Training)
PostHoc.Training2 <- emmeans(Hall2.Training, ~ Day,
                             adjust="bonferroni")
pairs(PostHoc.Training2)

PostHoc.Training3 <- emmeans(Hall2.Training, ~ Group*Day, adjust="bonferroni")
pairs(PostHoc.Training3)
cld(PostHoc.Training)

First.Hall.Training$Group <- factor(First.Hall.Training$Group, levels = c("Control", "Moderate warming", "High warming"))

ggboxplot(First.Hall.Training, x = "Day", y = "HallNumber",
          add = "jitter", shape = "Group",
          color = "black", fill="Group",
          xlab="Associative learning days",
          title="(a) ",
          ylab = ("Hall number"),
          palette = c("#00AFBB", "#E7B800", "#FC4E07")) +
  theme_bw() 



######################################## Fig 4b
#Fig 4b - learning and memory test

First.Hall.LM <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 5)

attach(First.Hall.LM)

Hall1 <- lmer(Number_first_hall ~ Group * Test + (1 | ID),
              data = First.Hall.LM)

Hall2 <- lmer(Number_first_hall ~ Group + Test + (1 | ID),
              data = First.Hall.LM)

Hall3 <- lmer(Latency ~ Number_first_hall + Group + Test + (1 | ID),
              data = First.Hall.LM)

AICctab(Hall1, Hall2, Hall3)
plotresid(Hall2, shapiro = T)
summary(Hall2)
Anova(Hall2, test.statistic="Chisq", type="III")
#poshoc 
phNH <- emmeans(Hall2, ~ Group, adjust="bonferroni")
pairs(phNH)
phNH2 <- emmeans(Hall2, ~ Test, adjust="bonferroni")
pairs(phNH2)
phNH3 <- emmeans(Hall2, ~ Test*Group, adjust="bonferroni")
pairs(phNH3)

First.Hall.LM$Group <- factor(First.Hall.LM$Group, levels = c("Control", "Moderate warming", "High warming"))

ggboxplot(First.Hall.LM, x = "Test", y = "Number_first_hall",
          add = "jitter", shape = "Group",
          color = "black", fill="Group", xlab="",
          title="(b) learning and memory test",
          ylab = ("Hall number"),
          palette = c("#00AFBB", "#E7B800", "#FC4E07")) +
  theme_bw()