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)

##################################################################
#Locomotor parameters 
Locomotor <- read_excel("raw data MS Impacts of ocean warming.xlsx", sheet = 6)

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

Locomotor$Test <- factor(Locomotor$Test, levels = c("Learning", "Memory"))

attach(Locomotor)

#Mean speed
descdist(`mean speed`, discrete=FALSE) #Gaussian
plot(density(`mean speed`))
hist(`mean speed`)

m1 <- glmer(`mean speed` ~ Group * Test + (1 | ID),
            data = Locomotor,
            family=Gamma(link = "log"))

plotresid(m1, shapiro = T)

summary(m1)
Anova(m1, test.statistic="Chisq", type="III")
em1 <- emmeans(m1, ~ Group)
pairs(em1)
cld(em1)

em2 <- emmeans(m1, ~ Day)
pairs(em2)
cld(em2)

em3 <- emmeans(m1, ~ Group * Test)
pairs(em3)
cld(em3)


Fig5a <- ggboxplot(Locomotor, x = "Group", y = "mean speed",
          ylab = "Average speed (cm/s)", xlab = "",
          add = "dotplot", shape = "Group",
          color = "black", ylim = c(0,15),
          fill = "Test", title="(a) ",
          palette = c('#8491A3','#0d7692')) +
  theme_bw() 
Fig5a



#maximum speed
descdist(`maximum speed`, discrete=FALSE) #Gaussian
plot(density(`maximum speed`))
hist(`mean speed`)

m2 <- glmer(`maximum speed` ~ Group * Test + (1 | ID),
            data = Locomotor,
            family=Gamma(link = "log"))

plotresid(m2, shapiro = T)

summary(m2)
Anova(m2, test.statistic="Chisq", type="III")


Fig5b <- ggboxplot(Locomotor, x = "Group", y = "maximum speed",
          ylab = "Maximum speed (cm/s)", xlab = "",
          add = "dotplot", shape = "Group",
          color = "black", 
          fill = "Test", title="(b) ",
          palette = c('#8491A3','#0d7692')) +
  theme_bw() 
Fig5b

#total distance
descdist(`total distance`, discrete=FALSE) #Gaussian

m3 <- glmer(`total distance` ~ Group * Test + (1 | ID),
            data = Locomotor,
            family=Gamma(link = "log"))

plotresid(m3, shapiro = T)

summary(m3)
Anova(m3, test.statistic="Chisq", type="III")

Fig5c <- ggboxplot(Locomotor, x = "Group", y = "total distance",
          ylab = "Distance traveled (cm)", xlab = "", 
          add = "dotplot", shape = "Group",
          color = "black", 
          fill = "Test", title="(c) ",
          palette = c('#8491A3','#0d7692')) +
  theme_bw()
Fig5c

#total time stopped
descdist(`total time stoped`, discrete=FALSE) #Gaussian

m4 <- glmer(`total time stoped` ~ Group * Test + (1 | ID),
            data = Locomotor,
            family=Gamma(link = "log"))

plotresid(m4, shapiro = T)

summary(m4)
Anova(m4, test.statistic="Chisq", type="III")
e4 <- emmeans(m4, ~ Group)
pairs(e4)
cld(e4)

em4 <- emmeans(m4, ~ Test)
pairs(em4)
cld(em4)

emt4 <- emmeans(m4, ~ Group * Test)
pairs(emt4)
cld(emt4)

Fig5d <- ggboxplot(Locomotor, x = "Group", y = "total time stoped",
          ylab = "Time stopped (s)", xlab = "", 
          add = "dotplot", shape = "Group",
          color = "black", 
          fill = "Test", title="(d) ",
          palette = c('#8491A3','#0d7692')) +
  theme_bw()
Fig5d

ggarrange(Fig5a, Fig5b, Fig5c, Fig5d, ncol = 2, nrow = 2,
          common.legend = T, legend = "bottom")