rm(list=ls())
library(tidyverse)
library(ggpubr)
library(knitr)
library(rstatix)
source("fun.R")

# species vectors for plotting:
cols <- c("AI" = "#FDE725FF", "KO" = "#5DC863FF","AC" = "#21908CFF",  "AV" = "#3B528BFF", "SA" = "#440154FF")
speciesnames <- c("AI" = "Acanthus ilicifolius", "KO" = "Kandelia obovata", "AC" = "Aegiceras corniculatum", "AV" = "Avicennia marina", "SA" = "Sonneratia apetala")

# load data
df.drag <- read.csv("dfdrag.csv")
df.drag$test<-factor(df.drag$test, levels = c("PC", "PW", "CW15"))
df.drag$species <-  factor(df.drag$species, levels = c("AI", "KO", "AC", "AV","SA"))

# make facetwrap labels
leaf.labs <- c("With leaves", "Without leaves")
names(leaf.labs) <- c("with", "without")
test.labs <- c("Current", "Wave", "Current + Wave")
names(test.labs) <- c("PC", "PW", "CW15")

### --- Drag force per surface area -----
# get regression
fit.c0_wl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "with", test == "PC"))
fit.w0_wl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "with", test == "PW"))
fit.cw_wl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "with", test == "CW15"))
fit.c0_nl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "without", test == "PC"))
fit.w0_nl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "without", test == "PW"))
fit.cw_nl <- lm(meanDrag ~ A_m2, filter(df.drag, leaves == "without", test == "CW15"))

# define cutoff values so lines do not extend beyond the data
sim_leafsize_wl <- seq(0.0149, 0.0652, 0.005)
sim_leafsize_nl <- seq(0.0037, 0.0228, 0.005)

# make regression lines
pred.c0_wl <- predict(fit.c0_wl, list(A_m2 = sim_leafsize_wl))
pred.w0_wl <- predict(fit.w0_wl, list(A_m2 = sim_leafsize_wl))
pred.cw_wl <- predict(fit.cw_wl, list(A_m2 = sim_leafsize_wl))
pred.c0_nl <- predict(fit.c0_nl, list(A_m2 = sim_leafsize_nl))
pred.w0_nl <- predict(fit.w0_nl, list(A_m2 = sim_leafsize_nl))
pred.cw_nl <- predict(fit.cw_nl, list(A_m2 = sim_leafsize_nl))

# make dataframe of regression lines
df.pred_drag <- data.frame(rbind(cbind(rep("PC", length(sim_leafsize_wl)),   rep("with", length(sim_leafsize_wl)), sim_leafsize_wl, pred.c0_wl),
                                 cbind(rep("PW", length(sim_leafsize_wl)),   rep("with", length(sim_leafsize_wl)), sim_leafsize_wl, pred.w0_wl),
                                 cbind(rep("CW15", length(sim_leafsize_wl)), rep("with", length(sim_leafsize_wl)), sim_leafsize_wl, pred.cw_wl),
                                 cbind(rep("PC", length(sim_leafsize_nl)),   rep("without", length(sim_leafsize_nl)), sim_leafsize_nl, pred.c0_nl),
                                 cbind(rep("PW", length(sim_leafsize_nl)),   rep("without", length(sim_leafsize_nl)), sim_leafsize_nl, pred.w0_nl),
                                 cbind(rep("CW15", length(sim_leafsize_nl)), rep("without", length(sim_leafsize_nl)), sim_leafsize_nl, pred.cw_nl)))
# adjust dataframe
colnames(df.pred_drag) <- c("test", "leaves", "A_m2", "drag")
df.pred_drag$A_m2 <- as.numeric(as.character(df.pred_drag$A_m2))
df.pred_drag$drag  <- as.numeric(as.character(df.pred_drag$drag))
df.drag$test<-factor(df.drag$test, levels = c("PC", "PW", "CW15"))
df.pred_drag$test<-factor(df.pred_drag$test, levels = c("PC", "PW", "CW15"))

# make plot
p.drag <-
  ggplot(filter(df.drag, !is.na(test), !is.na(leaves)),
         aes(x = A_m2, y = meanDrag, group = testname)) +
  geom_point(aes(colour = species)) +
  geom_line(data = df.pred_drag,
            aes(x = A_m2, y = drag, group = test)) +
  facet_grid(leaves~test, 
             labeller = labeller(leaves = leaf.labs, test = test.labs)) +
  scale_colour_manual(values = cols, labels = speciesnames) +
  labs(x = expression(paste("Surface area ",  A[proj], " (", m^2, ")")), 
       y = expression(paste("Maximum drag force ", F[D], " (N)")), 
       colour = "Species", title = "a") +
  theme_classic() +
  guides(colour = guide_legend(nrow = 2, byrow = FALSE))

### --- Cd per species -----
p.Cd_withleaves <-
  ggplot(filter(df.drag, !is.na(test), !is.na(leaves), leaves == "with"),
         aes(x = species, y = meanCd)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.2, aes(fill = species)) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.7) +
  scale_fill_manual(values = cols, labels = speciesnames)  +
  labs(x = "",
       y = expression(paste("Drag coefficient ", C[d])),
       fill = "Species", title = "b")+
  facet_wrap(~test,
             labeller = labeller(leaves = leaf.labs, test = test.labs)) +
  theme_classic() +
  theme(legend.position="bottom") 

p.Cd_withoutleaves <-
  ggplot(filter(df.drag, !is.na(test), !is.na(leaves), leaves == "without"),
         aes(x = species, y = meanCd)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.2, aes(fill = species)) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.7) +
  scale_fill_manual(values = cols, labels = speciesnames)  +
  labs(x = "Species",
       y = expression(paste("Drag coefficient ", C[d])),
       fill = "Species")+
  facet_wrap(~test,
             labeller = labeller(leaves = leaf.labs, test = test.labs)) +
  theme_classic() +
  theme(legend.position="bottom") +
  theme(strip.text.x = element_blank())

### --- plot all data -----
ggarrange(p.drag,
          p.Cd_withleaves, 
          p.Cd_withoutleaves, 
          nrow = 3, common.legend = TRUE, legend = "bottom", heights = c(2,1,1))

### --- statistics ----
#get information on regressions
summary(fit.c0_wl)
summary(fit.w0_wl)
summary(fit.cw_wl)
summary(fit.c0_nl)
summary(fit.w0_nl)
summary(fit.cw_nl)

# species comparisons
# check normality
df.drag %>%
  group_by(species, test, leaves) %>%
  shapiro_test(meanCd)
# check if there are species differences
df.drag %>%
  group_by(test, leaves) %>%
  kruskal_test(meanCd ~ species)
# check which species differences
df.drag %>%
  filter(test == "PC", leaves == "without") %>%
  dunn_test(meanCd ~ species)

### --- Cd table -----
df.drag %>%
  group_by(test, leaves, species) %>%
  summarise(mean_Cd = mean(Cd),
            SD_Cd = sd(Cd),
            SE_Cd = SD_Cd/sqrt(length(Cd))) %>%
  kable()