rm(list=ls())
library(tidyverse)
library(ggpubr)
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.branch <- data.frame(read.csv("dfbranch.csv", stringsAsFactors = FALSE, sep = ',', header = TRUE, fileEncoding = "UTF-8-BOM", na.strings = "NA"))
df.branch$species <- factor(df.branch$species, levels = c("AI", "KO", "AC", "AV","SA"))
df.branch$site <- factor(df.branch$site, levels = c("PR1", "PR2", "PR3", "PR4", "PR5", "ZH1", "HL1"))

### --- Maximum load per branch diameter -----
# make regression line
mod.sizeforce<- lm(formula = Fmax_N ~ I(diam_cm^2), data = filter(df.branch, arm_vs_diam == "correct"))
diams <- seq(min(df.branch[df.branch$arm_vs_diam == "correct",]$diam_cm), max(df.branch[df.branch$arm_vs_diam == "correct",]$diam_cm), 0.1)
pred.force <- predict(mod.sizeforce, list(diam_cm = diams))
df.predforce <- data.frame(diam_cm = diams, Fmax_N = pred.force)
df.predforce$species <- NA

# make plot
p.force_diam <- ggplot(filter(df.branch, arm_vs_diam == "correct"), 
                              aes(x = diam_cm, y = Fmax_N, group = species)) + 
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  labs(x = expression(paste("Branch diameter ", symbol("\306"), " (cm)")),
       y = expression(paste("Maximum load ", F[max], " (N)")) , 
       colour = "Species",
       title = "a") +
  scale_x_continuous(breaks = seq(0, 2, 0.5), limits = c(0, 2)) +
  theme_classic() +
  geom_line(data = df.predforce, aes(x = diam_cm, y = Fmax_N), colour = "black") +
  theme(legend.position="bottom") +
  guides(colour = guide_legend(nrow = 2, byrow = FALSE))

### --- Force needed to bend branch, per branch diameter -----
# make regression line
mod.sizebend<- lm(formula = slope_N.mm ~ I(diam_cm^2), 
                  data = filter(df.branch, arm_vs_diam == "correct"))
diams <- seq(min(df.branch[df.branch$arm_vs_diam == "correct",]$diam_cm), max(df.branch[df.branch$arm_vs_diam == "correct",]$diam_cm), 0.1)
pred.bend <- predict(mod.sizebend, list(diam_cm = diams))
df.predbend <- data.frame(diam_cm = diams, slope_N.mm = pred.bend)
df.predbend$species <- NA

# make plot
p.bend_diam <- ggplot(filter(df.branch, arm_vs_diam == "correct"), 
                             aes(x = diam_cm, y = slope_N.mm, group = species)) + 
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  labs(x = expression(paste("Branch diameter ", symbol("\306"), " (cm)")),
       y = expression(paste("Force to bend F/x (N ", mm^-1, ")")), 
       colour = "Species",
       title = "b") +
  scale_x_continuous(breaks = seq(0, 2, 0.5), limits = c(0, 2)) +
  theme_classic() +
  geom_line(data = df.predbend, aes(x = diam_cm, y = slope_N.mm), colour = "black") 

### --- Branch strength per species for branch diameters 1.25-1.75 cm -----
# prepare data
df.subbranch <- 
  df.branch %>%
  filter(arm_vs_diam == "correct",diam_cm > 1.25, diam_cm < 1.75) %>%
  group_by(ID) %>%
  mutate(mean_diam_cm = mean(diam_cm, na.rm = TRUE),
         mean_Fmax_N  = mean(Fmax_N, na.rm = TRUE)) %>%
  distinct(ID, .keep_all = TRUE)

# make plot
p.strength_species <- 
  ggplot(data = df.subbranch, 
         aes(x = species, y = Fmax_N)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.6) +
  geom_violin(aes(fill = species, alpha = 0.8)) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  scale_fill_manual(values = cols,
                    labels = speciesnames) +
  labs(x = "Species", y = "Max. load 1.25-1.75 cm (N)", fill = "Species", title = "c") +
  theme_classic() 

### --- Modulus of Rupture per species -----
# make plot
p.MOR_species <- 
  ggplot(df.branch, 
         aes(x = species, y = meanMOR_MPa)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.6) +
  geom_violin(aes(fill = species, alpha = 0.8)) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  scale_fill_manual(values = cols,
                    labels = speciesnames) +
  labs(x = "Species", 
       y = expression(paste("MOR (N ", mm^-2, ")")),
       fill = "Species",
       title = "d") +
  theme_classic() +
  scale_alpha(guide = 'none') +
  scale_y_continuous(breaks = seq(0, 70, 10), limits = c(0, 65)) 

### --- Modulus of Elasticity per species -----
# make plot
p.MOE_species <- 
  ggplot(df.branch, 
         aes(x = species, y = meanMOE_MPa)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.3) +
  geom_violin(aes(fill = species, alpha = 0.8)) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  scale_fill_manual(values = cols,
                    labels = speciesnames) +
  labs(x = "Species", 
       y = expression(paste("MOE (N ", mm^-2, ")")),
       fill = "Species",
       title = "e") +
  theme_classic() +
  scale_alpha(guide = 'none') +
  scale_y_continuous(breaks = seq(0, 6000, 1000), limits = c(0, 5300))

### --- MOR vs MOE -----
# make regression line
mod.mormoe<- lm(meanMOE_MPa ~ meanMOR_MPa, df.branch)
mors <- seq(min(df.branch$meanMOR_MPa), max(df.branch$meanMOR_MPa), 0.1)
pred.moe <- predict(mod.mormoe, list(meanMOR_MPa = mors))
df.predmoe <- data.frame(meanMOR_MPa = mors, meanMOE_MPa = pred.moe)
df.predmoe$species <- NA

# make plot
p.MOR_MOE <- 
  ggplot(data = df.branch, aes(x = meanMOR_MPa, y = meanMOE_MPa)) + 
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  scale_x_continuous(breaks = seq(0, 70, 10), limits = c(0, 62)) +
  scale_y_continuous(breaks = seq(0, 6000, 1000), limits = c(0, 5300)) +
  labs(x = expression(paste("MOR (N ", mm^-2, ")")),
       y = expression(paste("MOE (N ", mm^-2, ")")), 
       colour = "Species", title = "f") +
  theme_classic() +
  geom_line(data = df.predmoe, aes(x = meanMOR_MPa, y = meanMOE_MPa), colour = "black") 

### --- plot data -----
ggarrange(p.force_diam, p.bend_diam,
          p.strength_species, p.MOR_species,
          p.MOE_species, p.MOR_MOE,
          ncol = 2, nrow = 3,
          common.legend = TRUE, legend = "bottom")

### --- statistics ----
#get information on regressions
summary(mod.sizeforce)
summary(mod.sizebend)
summary(mod.mormoe)

# species comparisons
# check normality
df.subbranch %>%
  group_by(species) %>%
  shapiro_test(mean_Fmax_N)
df.branch %>%
  group_by(species) %>%
  shapiro_test(meanMOR_MPa)
df.branch %>%
  group_by(species) %>%
  shapiro_test(meanMOE_MPa)

# check if there are species differences
df.subbranch %>%
  ungroup() %>%
  kruskal_test(mean_Fmax_N  ~ species) 
df.branch %>%
  kruskal_test(meanMOR_MPa ~ species)
df.branch %>%
  kruskal_test(meanMOE_MPa ~ species)

# check which species differences
df.subbranch %>%
  ungroup() %>%
  dunn_test(mean_Fmax_N ~ species) %>%
  add_significance()
df.branch %>%
  dunn_test(meanMOR_MPa ~ species)
df.branch %>%
  dunn_test(meanMOE_MPa ~ species)