##### Morphology PeerJ

library(DescTools)
library(dplyr)
library(lmtest)
library(rstatix)
library(ggpubr)

### Kruskal-wallis for each site

leafl <- read.csv("leaflen.csv", sep=",")
leafw <- read.csv("leafwid.csv", sep=",")

str(leafl)
leafl$site
leafl$lenght
leafw$site
leafw$width

leafl$site<-factor(leafl$site)
leafw$site<-factor(leafw$site)

plot (lenght~site, leafl, ylab="Leaf lenght (cm)", xlab="Sites")
plot (width~site, leafw,  ylab="Leaf width (mm)", xlab="Sites")

ggboxplot(leafl, x = "site", y = "lenght", 
          color = "site", palette = c("black", "purple", "dark green", "orange", "brown"),
          order = c("PS", "MAC", "SETU", "PM", "LC"),
          ylab = "Leaf lenght  (cm)", xlab = "Sites")

ggboxplot(leafw, x = "site", y = "width", 
          color = "site", palette = c("black", "purple", "dark green", "orange", "brown"),
          order = c("PS", "MAC", "SETU", "PM", "LC"),
          ylab = "Leaf width  (mm)", xlab = "Sites")

###Leaf lenght

### Leaf lengh
kruskal.test(lenght ~ site, leafl)

##Dunn post hoc
DunnTest(lenght ~ site, leafl, method = "bonferroni")

###Leaf width
kruskal.test(width ~ site, leafw)

##Dunn post hoc
DunnTest(width ~ site, leafw, method ="bonferroni")

## Effect size

leafl %>% kruskal_effsize(lenght ~ site)

leafw %>% kruskal_effsize(width ~ site)

###### Comparisons by leaf tip shape

lenshape <- read.csv("shapelen.csv", header=T, dec=".", sep=",")
widshape <- read.csv("shapewid.csv", header=T, dec=".", sep=",")
lenshape$leaf
lenshape$lenght
widshape$leaf
widshape$width

lenshape$leaf<-factor(lenshape$leaf)
is.factor(lenshape$leaf)
plot (lenght~leaf, lenshape)

widshape$leaf<-factor(widshape$leaf)
is.factor(widshape$leaf)
plot (width~leaf, widshape)

ggboxplot(lenshape, x = "leaf", y = "lenght", 
          color = "leaf", palette = c("black", "orange"),
          order = c("HB", "HW"),
          ylab = "Leaf lenght  (cm)", xlab = "Leaf-tip")

ggboxplot(widshape, x = "leaf", y = "width", 
          color = "leaf", palette = c("black", "orange"),
          order = c("HB", "HW"),
          ylab = "Leaf width  (mm)", xlab = "Leaf-tip")

### Leaf lenght
#### Mann-whitney u leaf lenght

wilcox.test(lenght ~ leaf, lenshape, exact = FALSE)

### Leaf width
#### Mann-whitney u leaf lenght

wilcox.test(width ~ leaf, widshape, exact = FALSE)