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

### --- Leaf surface area per species -----
# make plot
p.area_species <- 
  ggplot(df.leaf, aes(x = species, y = area_cm2)) +
  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(A[leaf], " (", cm^2, ")")), fill = "Species",
       title = "a") +
  theme_classic() +
  scale_alpha(guide = 'none') +
  theme(legend.position="bottom") +
  guides(fill = guide_legend(nrow = 2, byrow = FALSE))

### --- Pulling force per species -----
# make plot
p.force_species <- 
  ggplot(df.leaf, aes(x = species, y = Fmax_N)) +
  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(F[pull], " (N)")), fill = "Species",
       title = "b") +
  theme_classic() +
  scale_alpha(guide = 'none') +
  scale_y_continuous(breaks = seq(0, 27, 5), limits = c(0, 25)) 

### --- Leaf surface area vs pulling force
# make regression line
mod.area_pull <- lm(Fmax_N ~ area_cm2, df.leaf)
areas <- seq(min(df.leaf$area_cm2, na.rm = TRUE), max(df.leaf$area_cm2, na.rm = TRUE), 0.1)
pred.pull <- predict(mod.area_pull, list(area_cm2 = areas))
df.predpull <- data.frame(area_cm2 = areas, Fmax_N = pred.pull)
df.predpull$species <- NA

# make plot
summary(mod.area_pull)
p.area_force <- 
  ggplot(df.leaf, aes(x = area_cm2, y = Fmax_N, group = species)) +
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  labs(x = expression(paste(A[leaf], " (", cm^2, ")")), 
       y = expression(paste(F[pull], " (N)")), 
       colour = "Species", title = "c") +
  theme_classic() +
  scale_x_continuous(breaks = seq(0, 70, 20), limits = c(0, 70)) +
  scale_y_continuous(breaks = seq(0, 25, 5), limits = c(0, 25)) +
  geom_line(data = df.predpull, aes(x = area_cm2, y = Fmax_N), colour = "black") 

### --- Petiole diameter vs pulling force
# make regression line
mod.diam_pull <- lm(Fmax_N ~ petioleDiam_cm, df.leaf)
diams <- seq(min(df.leaf$petioleDiam_cm, na.rm = TRUE), max(df.leaf$petioleDiam_cm, na.rm = TRUE), 0.01)
pred.diampull <- predict(mod.diam_pull, list(petioleDiam_cm = diams))
df.preddiampull <- data.frame(petioleDiam_cm = diams, Fmax_N = pred.diampull)
df.preddiampull$species <- NA

# make plot
p.diam_force <- 
  df.leaf %>%
  filter(!is.na(petioleDiam_cm)) %>%
  ggplot(aes(x = petioleDiam_cm, y = Fmax_N, group = species)) +
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  theme_classic() +
  labs(x = "Petiole diameter (cm)", 
       y = expression(paste(F[pull], " (N)")),
       colour = "Species",
       title = "d") +
  scale_x_continuous(limits = c(0, 0.5)) +
  scale_y_continuous(limits = c(0, 25)) +
  geom_line(data = df.preddiampull, aes(x = petioleDiam_cm, y = Fmax_N), colour = "black") +
  guides(colour = guide_legend(nrow = 2, byrow = FALSE))

### --- Modulus of Leaf loss (area-based) per species -----
# make plot
p.MOLA_species <- 
  ggplot(df.leaf, aes(x = species, y = (Fmax_N/(area_cm2)))) +
  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(MOL[A], " (N ", cm^-2,")")), fill = "Species",
       title = "e") +
  theme_classic() +
  scale_alpha(guide = 'none') 

### --- Leaf Mass per Area per species -----
# make plot
p.LMA_species <-
  ggplot(df.leaf, aes(x = species, y = lma_gcm2)) +
  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("LMA (", g, " ", cm^-2, ")")), fill = "Species",
       title = "f") +
  theme_classic() +
  scale_alpha(guide = 'none') 

### --- Leaf Mass per Area per pulling force
# make plot
summary(lm(Fmax_N ~ lma_gcm2, df.leaf))
p.LMA_force <- 
  ggplot(df.leaf, aes(x = lma_gcm2, y = Fmax_N, group = species)) +
  geom_point(aes(colour = species)) +
  scale_colour_manual(values = cols,
                      labels = speciesnames) +
  labs(x = expression(paste("LMA (", g, " ", cm^-2, ")")),
       y = expression(paste(F[pull], " (N)")), 
       colour = "Species", title = "g") +
  theme_classic() +
  scale_x_continuous(breaks = seq(0, 0.03, 0.01), limits = c(0, 0.03)) +
  scale_y_continuous(breaks = seq(0, 25, 5), limits = c(0, 25)) 

### --- Modulus of Leaf loss (LMA-based) per species -----
# make plot
p.MOLM_species <- 
  ggplot(df.leaf,
         aes(x = species, y = Fmax_N/lma_gcm2)) +
  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(MOL[M], " (N  ", g^-1, " ", cm^2, ")")), 
       fill = "Species",
       title = "h") +
  theme_classic() 

### --- plot data -----
ggarrange(p.area_species, p.force_species,
          p.area_force, p.diam_force,
          p.MOLA_species, p.LMA_species, 
          p.LMA_force, p.MOLM_species,
          ncol = 2, nrow = 4,
          common.legend = TRUE, legend = "bottom")

### --- statistics ----
summary(mod.area_pull)
summary(mod.diam_pull)

# species comparisons
# check normality
df.leaf %>%
  group_by(species) %>%
  shapiro_test(area_cm2)
df.leaf %>%
  group_by(species) %>%
  shapiro_test(Fmax_N)
df.leaf %>%
  group_by(species) %>%
  shapiro_test(lma_gcm2)
df.leaf %>%
  mutate(MOLA = Fmax_N/area_cm2) %>%
  group_by(species) %>%
  shapiro_test(MOLA)
df.leaf %>%
  mutate(MOLM = Fmax_N/lma_gcm2) %>%
  group_by(species) %>%
  shapiro_test(MOLM)

# check if there are species differences
df.leaf %>%
  kruskal_test(area_cm2 ~ species)
df.leaf %>%
  kruskal_test(Fmax_N ~ species)
df.leaf %>%
  kruskal_test(lma_gcm2 ~ species)
df.leaf %>%
  mutate(MOLA = Fmax_N/area_cm2) %>%
  kruskal_test(MOLA ~ species)
df.leaf %>%
  mutate(MOLM = Fmax_N/lma_gcm2) %>%
  kruskal_test(MOLM ~ species)

# check which species differences
df.leaf %>%
  dunn_test(area_cm2 ~ species)
df.leaf %>%
  dunn_test(Fmax_N ~ species)
df.leaf %>%
  dunn_test(lma_gcm2 ~ species)
df.leaf %>%
  mutate(MOLA = Fmax_N/area_cm2) %>%
  dunn_test(MOLA ~ species)
df.leaf %>%
  mutate(MOLM = Fmax_N/lma_gcm2) %>%
  dunn_test(MOLM ~ species)