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 2 - Target vs Halls - sheet 1
Maze.Halls <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 1)

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

attach(Maze.Halls)
descdist(Time, discrete=FALSE) 
plot(density(Time))
hist(Time)

TimeL <- sqrt(Time)
descdist(TimeL, discrete=FALSE) 
plot(density(TimeL))
hist(TimeL)

### Fig 2a
# Control
datControl <-  Maze.Halls %>%
  filter(Group == "Control") %>%
  mutate(Time_squared = Time^2)

c1 <- glmer(Time_squared ~ Test * `Maze arm` + (1 | ID),
            data = datControl,
            family=Gamma(link = "log"))

c2 <- glmer(Time_squared ~ Test + `Maze arm` + (1 | ID),
            data = datControl,
            family=Gamma(link = "log"))

AICctab(c1, c2)
plotresid(c1, shapiro = T)
plotresid(c2, shapiro = T)
summary(c2)
Anova(c2, test.statistic="Chisq", type="III")
em1 <- emmeans(c2, ~ `Maze arm`, adjust = "bonferroni")
pairs(em1)
cld(em1)


### Fig 2b
# Moderate
datModerate <-  Maze.Halls %>%
  filter(Group == "Moderate warming") %>%
  mutate(Time_squared = Time^2)

m1 <- glmer(Time_squared ~ Test * `Maze arm` + (1 | ID),
            data = datModerate,
            family=Gamma(link = "log"))

m2 <- glmer(Time_squared ~ Test + `Maze arm` + (1 | ID),
            data = datModerate,
            family=Gamma(link = "log"))

AICctab(m1, m2)
plotresid(m1, shapiro = T)
plotresid(m2, shapiro = T)
summary(m2)
Anova(m2, test.statistic="Chisq", type="III")
em2 <- emmeans(m2, ~ `Maze arm`, adjust = "bonferroni")
pairs(em2)
cld(em2)

#Fig 2c
# High
datHigh <-  Maze.Halls %>%
  filter(Group == "High warming") %>%
  mutate(Time_squared = Time^2)

h1 <- glmer(Time_squared ~ Test * `Maze arm` + (1 | ID),
            data = datHigh,
            family=Gamma(link = "log"))

h2 <- glmer(Time_squared ~ Test + `Maze arm` + (1 | ID),
            data = datHigh,
            family=Gamma(link = "log"))

AICctab(h1, h2)
plotresid(h1, shapiro = T)
plotresid(h2, shapiro = T)
summary(h2)
Anova(h2, test.statistic="Chisq", type="III")
em3 <- emmeans(h2, ~ `Maze arm`, adjust = "bonferroni")
pairs(em3)
cld(em3)

#####Fig2a
Fig2a <- ggboxplot(datControl, x ="Maze arm", y = "Time",
          ylab = "Time spent in different areas of the tank (%)",
          xlab = "", 
          add = "dotplot", shape = "Group",
          color = "black", ylim=c(0,100),
          fill = "Test", title="(a) Control ",
          palette = c('#8491A3','#0d7692')) +
  theme_bw()
Fig2a

#####Fig2b
Fig2b <- ggboxplot(datModerate, x ="Maze arm", y = "Time",
                   ylab = "Time spent in different areas of the tank (%)",
                   xlab = "", ylim=c(0,100),
                   add = "dotplot", shape = "Group",
                   color = "black", 
                   fill = "Test", title="(b) Moderate warming ",
                   palette = c('#8491A3','#0d7692')) +
  theme_bw()
Fig2b



#####Fig2c
Fig2c <- ggboxplot(datHigh, x ="Maze arm", y = "Time",
                   ylab = "Time spent in different areas of the tank (%)",
                   xlab = "", 
                   add = "dotplot", shape = "Group",
                   color = "black", ylim=c(0,100),
                   fill = "Test", title="(c) High warming",
                   palette = c('#8491A3','#0d7692')) +
  theme_bw()
Fig2c

ggarrange(Fig2a, Fig2b, Fig2c, ncol = 3, nrow = 1,
          common.legend = T, legend = "bottom")