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"))

df.leaf <- data.frame(read.csv("dfleaf.csv", stringsAsFactors = FALSE, sep = ',', header = TRUE, fileEncoding = "UTF-8-BOM", na.strings = "NA"))
df.leaf$species <- factor(df.leaf$species, levels = c("AI", "KO", "AC", "AV","SA"))
df.leaf$site <- factor(df.leaf$site, levels = c("PR1", "PR2", "PR3", "PR4", "PR5", "ZH1", "HL1"))
colnames(df.leaf) <- c("species", "site", "tree", "ID", "Fpull_N", "area_cm2", "dryweight_g", "lma_gcm2")

# prepare dataframes
df.branchleaf <- df.branch %>%
  filter(arm_vs_diam == "correct", diam_cm > 1.25, diam_cm < 1.75) %>%
  left_join(df.leaf[c("ID", "Fpull_N")], by = c("ID"))
df.branchleaf$species <- factor(df.branchleaf$species, levels = c("AI", "KO", "AC", "AV","SA"))

df.meanforce <- df.branchleaf %>%
  group_by(species) %>%
  summarise(muBranch = mean(na.omit(Fmax_N)), muLeaf = mean(na.omit(Fpull_N)))
df.meanforce$species <- factor(df.meanforce$species, levels = c("AI", "KO", "AC", "AV","SA"))

### --- plot branches max force versus leaves pulling force -----
# make plot
p.branch_leaf <-
  ggplot(filter(df.branchleaf, species %in% c("AI", "KO", "AC", "AV", "SA"), !is.na(species)), 
         aes(x = Fmax_N, y = Fpull_N, group = species)) +
  geom_point(data = df.meanforce, 
             aes(x = muBranch, y = muLeaf, colour = species), 
             size = 3, shape = 15) +
  geom_point(aes(colour = species), shape = 16, alpha = 0.7) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  labs(x = expression(paste("Branch maximum load ",  F[max] ,"(N)")), 
       y = expression(paste("Leaf detachment force ", F[pull] , "(N)")), 
       colour = "Species", title = "a") +
  theme_classic() +
  theme(legend.position="bottom") +
  guides(colour = guide_legend(nrow = 2, byrow = FALSE)) 

### --- plot strength ratio branches vs eaves -----
p.branch_leaf_ratio <-
  ggplot(filter(df.branchleaf, species %in% c("AI", "KO", "AC", "AV", "SA")),
         aes(x = species, y = Fpull_N/Fmax_N)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.45) +
  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) +
  theme_classic() +
  labs(x = "Species", 
       y = expression(paste(F[pull], " / ", F[max], " (N/N)")), 
       colour = "Species", title = "b") +
  scale_y_log10()

### --- plot data -----
ggarrange(p.branch_leaf, p.branch_leaf_ratio, common.legend = TRUE, legend = "bottom")

### --- statistics ----
# species comparisons
# check normality
df.branchleaf %>%
  mutate(ratio = Fpull_N/Fmax_N) %>%
  group_by(species) %>%
  shapiro_test(ratio)
# check if there are species differences
df.branchleaf %>%
  mutate(ratio = Fpull_N/Fmax_N) %>%
  kruskal_test(ratio ~ species)
# check which species differences
df.branchleaf %>%
  mutate(ratio = Fpull_N/Fmax_N) %>%
  dunn_test(ratio ~ species)