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

### --- Branch MOR per site per species ------
p.MOR_species <- 
  ggplot(df.branch, 
         aes(x = species, y = meanMOR_MPa)) + 
  geom_violin(aes(fill = species)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.3) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  stat_summary(aes(y = meanMOR_MPa, group = 1), fun = mean, colour = "red", geom = "line", group = 1) +
  facet_wrap(~site, nrow = 1)+
  scale_fill_manual(values = cols, labels = speciesnames) +
  labs(x = "", y = "MOR (MPa)", fill = "Species", title = "a") +
  theme_classic() +
  theme(legend.position="bottom") +
  guides(fill = guide_legend(nrow = 2, byrow = FALSE))

### --- Leaf Fpull per site per species ------
p.Fpull_species <- 
  ggplot(df.leaf, 
         aes(x = species, y = Fpull_N)) + 
  geom_violin(aes(fill = species)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.3) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  stat_summary(aes(y = Fpull_N, group = 1), fun = mean, colour = "red", geom = "line", group = 1) +
  facet_wrap(~site, nrow = 1)+
  scale_fill_manual(values = cols, labels = speciesnames) +
  labs(x = "Species", y = "Fpull (N)", fill = "Species") +
  theme_classic() +
  theme(strip.text.x = element_blank())

### --- Branch MOR per species per site ------
p.MOR_site <-
  ggplot(df.branch, 
       aes(x = site, y = meanMOR_MPa)) + 
  geom_violin(aes(fill = species)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.3) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  stat_summary(aes(y = meanMOR_MPa, group = 1), fun = mean, colour = "red", geom = "line", group = 1) +
  facet_wrap(~species, nrow = 1)+
  scale_fill_manual(values = cols, labels = speciesnames) +
  labs(x = "", y = "MOR (MPa)", fill = "Species", title = "b") +
  theme_classic() +
  theme(legend.position="bottom") +
  guides(fill = guide_legend(nrow = 2, byrow = FALSE))

### --- Leaf Fpull per species per site ------
p.Fpull_site <- 
  ggplot(df.leaf, 
       aes(x = site, y = Fpull_N)) + 
  geom_violin(aes(fill = species)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 0.3) +
  stat_summary(fun.data = data_summary, geom = "crossbar", width = 0.1, color = "red", fill = "white", alpha = 0.8) +
  stat_summary(aes(y = Fpull_N, group = 1), fun = mean, colour = "red", geom = "line", group = 1) +
  facet_wrap(~species, nrow = 1)+
  scale_fill_manual(values = cols, labels = speciesnames) +
  labs(x = "Site (ordered by salinity)", y = "Fpull (N)", fill = "Species") +
  theme_classic() +
  theme(strip.text.x = element_blank())

# plot data
ggarrange(p.MOR_species, p.Fpull_species, 
          p.MOR_site, p.Fpull_site, 
          common.legend = TRUE, nrow = 4)

### --- statistics -----
# per species, per site:
st <- "PR1" #fill in site of interest
# check normality
df.branch %>%
  filter(site == st) %>%
  group_by(species) %>%
  shapiro_test(meanMOR_MPa)
df.leaf %>%
  filter(site == st) %>%
  group_by(species) %>%
  shapiro_test(Fpull_N)
# check if there are differences
df.branch %>%
  filter(site == st) %>%
  kruskal_test(meanMOR_MPa ~ species)
df.leaf %>%
  filter(site == st) %>%
  kruskal_test(Fpull_N ~ species)
# check which differences
df.branch %>%
  filter(site == st) %>%
  dunn_test(meanMOR_MPa ~ species)
df.leaf %>%
  filter(site == st) %>%
  dunn_test(Fpull_N ~ species)

# per site, per species
sp <- "AI" #fill in species of interest
# check normality
df.branch %>%
  filter(species == sp) %>%
  group_by(site) %>%
  shapiro_test(meanMOR_MPa)
df.leaf %>%
  filter(species == sp) %>%
  group_by(site) %>%
  shapiro_test(Fpull_N)
# check if there are differences
df.branch %>%
  filter(species == sp) %>%
  kruskal_test(meanMOR_MPa ~ site)
df.leaf %>%
  filter(species == sp) %>%
  kruskal_test(Fpull_N ~ site)
# check which differences
df.branch %>%
  filter(species == sp) %>%
  dunn_test(meanMOR_MPa ~ site)
df.leaf %>%
  filter(species == sp) %>%
  dunn_test(Fpull_N ~ site)