# R code
# All code used is given for the submission at the top

#### Packages ####
library(carData)
library(car)
library(effects)
library(MASS)
library(ggplot2)
library(nortest)
library(Matrix)
library(lme4)
library(devtools)
library(rstatix)
library(sjPlot)
library(tidyverse)
library(GGally)
library(readr)
library(lsmeans)
library(emmeans)
library(sandwich)
library(glmmTMB)
library(heplots)
library(dunn.test)
library(gamlss)
library(DHARMa)
library(labdsv)
library(dplyr)
library(MASS)

LGSH <- read_csv("/Volumes/Kate_/Reefscapers/Spawning Research /LGSH_Final/Items to submit/Resubmission1/Submit/Supplementary File 1.csv")

#### Define Variables ####
LGSH$Spawn_Date <- as.Date(as.character(LGSH$Spawn_Date),format = "%Y%m%d")
LGSH$Pale <- as.Date(as.character(LGSH$Pale),format = "%Y%m%d")
LGSH$Pigmented <- as.Date(as.character(LGSH$Pigmented),format = "%Y%m%d")
LGSH$Eggs_Located <- as.Date(as.character(LGSH$Eggs_Located),format = "%Y%m%d")
LGSH$Reef_type <- as.factor(LGSH$Reef_type)
LGSH$Atoll <- as.factor(LGSH$Atoll)
LGSH$Year <- as.factor(LGSH$Year)
LGSH$Species <- as.factor(LGSH$Species)
LGSH$Genus <- as.factor(LGSH$Genus)
LGSH$Spawn_Month <- as.factor(LGSH$Spawn_Month)
LGSH$Days_after_full_moon <- as.numeric(LGSH$Days_after_full_moon)
LGSH$WindSpeed_DailyAV <- as.numeric(LGSH$WindSpeed_DailyAV)
LGSH$Total_Daily_Precipitation_mm <- as.numeric(LGSH$Total_Daily_Precipitation_mm)
LGSH$AverageDaily_SST <- as.numeric(LGSH$AverageDaily_SST)
LGSH$Tide_Depth_m <- as.numeric(LGSH$Tide_Depth_m)
LGSH$Moon_Illumination <- as.numeric(LGSH$Moon_Illumination)
LGSH$Proximity_to_FM <- as.integer(LGSH$Proximity_to_FM)
LGSH$Spawning_Time_After_Sunset_Minutes <- as.numeric(LGSH$Spawning_Time_After_Sunset_Minutes)
LGSH$Bundling_Time_After_Sunset_Minutes <- as.numeric(LGSH$Bundling_Time_After_Sunset_Minutes)

#### NAs ####
na.omit(LGSH) 

#### Subset by Species ####
# not all have observations from both atolls & some have very small sample sizes 
asperaLGSH <- subset(LGSH, Species=="aspera")
cythereaLGSH <- subset(LGSH, Species=="cytherea")
digitiferaLGSH <- subset(LGSH, Species=="digitifera")
gemmiferaLGSH <- subset(LGSH, Species=="gemmifera")
hempriciiLGSH <- subset(LGSH, Species=="hempricii")
humilisLGSH <- subset(LGSH, Species=="humilis")
hyacinthusLGSH <- subset(LGSH, Species=="hyacinthus")
milleporaLGSH <- subset(LGSH, Species=="millepora")
muricataLGSH <- subset(LGSH, Species=="muricata")
nasutaLGSH <- subset(LGSH, Species=="nasuta")
plantaLGSH <- subset(LGSH, Species=="plantaginea")
samoensisLGSH <- subset(LGSH, Species=="samoensis")
secaleLGSH <- subset(LGSH, Species=="secale")
squarrosaLGSH <- subset(LGSH, Species=="squarrosa")
tenuisLGSH <- subset(LGSH, Species=="tenuis")
validaLGSH <- subset(LGSH, Species=="valida")


#### testing response ###
# for all data 
nTL = length(LGSH$Proximity_to_FM)
shapiro.test(LGSH$Proximity_to_FM) # normal
leveneTest(Proximity_to_FM ~ Atoll, data = LGSH) # uncommon variance

nTL = length(LGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(LGSH$Spawning_Time_After_Sunset_Minutes) # normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = LGSH) # uncommon variance


# for all species
# prox FM
nTL = length(asperaLGSH$Proximity_to_FM)
shapiro.test(asperaLGSH$Proximity_to_FM) # normal
leveneTest(Proximity_to_FM ~ Atoll, data = asperaLGSH) # uncommon variance

nTL = length(cythereaLGSH$Proximity_to_FM)
shapiro.test(cythereaLGSH$Proximity_to_FM) # not normal

nTL = length(digitiferaLGSH$Proximity_to_FM)
shapiro.test(digitiferaLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = digitiferaLGSH) # common variance

nTL = length(gemmiferaLGSH$Proximity_to_FM)
shapiro.test(gemmiferaLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = gemmiferaLGSH) # common variance

nTL = length(hempriciiLGSH$Proximity_to_FM)
shapiro.test(hempriciiLGSH$Proximity_to_FM) # not normal

nTL = length(humilisLGSH$Proximity_to_FM)
shapiro.test(humilisLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = humilisLGSH) # common variance

nTL = length(hyacinthusLGSH$Proximity_to_FM)
shapiro.test(hyacinthusLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = hyacinthusLGSH) # common variance

nTL = length(milleporaLGSH$Proximity_to_FM)
shapiro.test(milleporaLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = milleporaLGSH) # uncommon variance

nTL = length(muricataLGSH$Proximity_to_FM)
shapiro.test(muricataLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = muricataLGSH) # common variance

nTL = length(nasutaLGSH$Proximity_to_FM)
shapiro.test(nasutaLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = nasutaLGSH) # common variance

nTL = length(plantaLGSH$Proximity_to_FM)
shapiro.test(plantaLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = plantaLGSH) # uncommon variance

nTL = length(samoensisLGSH$Proximity_to_FM)
shapiro.test(samoensisLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = samoensisLGSH) # uncommon variance

nTL = length(secaleLGSH$Proximity_to_FM)
shapiro.test(secaleLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = secaleLGSH)  # common variance

nTL = length(squarrosaLGSH$Proximity_to_FM)
shapiro.test(squarrosaLGSH$Proximity_to_FM) # normal
leveneTest(Proximity_to_FM ~ Atoll, data = squarrosaLGSH) # common variance

nTL = length(tenuisLGSH$Proximity_to_FM)
shapiro.test(tenuisLGSH$Proximity_to_FM) # not normal
leveneTest(Proximity_to_FM ~ Atoll, data = tenuisLGSH) # common variance

nTL = length(validaLGSH$Proximity_to_FM)
shapiro.test(validaLGSH$Proximity_to_FM) # normal
leveneTest(Proximity_to_FM ~ Atoll, data = validaLGSH) # common variance


# spawn time after SS
nTL = length(asperaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(asperaLGSH$Spawning_Time_After_Sunset_Minutes) # normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = asperaLGSH) # uncommon variance

nTL = length(cythereaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(cythereaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal

nTL = length(digitiferaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(digitiferaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = digitiferaLGSH) # common variance

nTL = length(gemmiferaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(gemmiferaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = gemmiferaLGSH) # common variance

nTL = length(hempriciiLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(hempriciiLGSH$Spawning_Time_After_Sunset_Minutes) # not normal

nTL = length(humilisLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(humilisLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = humilisLGSH) # common variance

nTL = length(hyacinthusLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(hyacinthusLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = hyacinthusLGSH) # common variance

nTL = length(milleporaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(milleporaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = milleporaLGSH) # uncommon variance

nTL = length(muricataLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(muricataLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = muricataLGSH) # common variance

nTL = length(nasutaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(nasutaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = nasutaLGSH) # common variance

nTL = length(plantaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(plantaLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = plantaLGSH) # uncommon variance

nTL = length(samoensisLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(samoensisLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = samoensisLGSH) # common variance

nTL = length(secaleLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(secaleLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = secaleLGSH)  # common variance

nTL = length(squarrosaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(squarrosaLGSH$Spawning_Time_After_Sunset_Minutes) # normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = squarrosaLGSH) # common variance

nTL = length(tenuisLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(tenuisLGSH$Spawning_Time_After_Sunset_Minutes) # not normal
leveneTest(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = tenuisLGSH) # common variance

nTL = length(validaLGSH$Spawning_Time_After_Sunset_Minutes)
shapiro.test(validaLGSH$Spawning_Time_After_Sunset_Minutes) # sample size too small


#### Kruskal Wallis tests ####
# don't assume normality or common variance 
library(rstatix)
# kruskal.test(response ~ explan, data = LGSH)
# kruskal_effsize(response ~ explan, data = LGSH) # gives eta^2
# f = sqrt(eta/(1-eta)) ; need to do this if using Cohen classification

kruskal.test(Proximity_to_FM ~ Atoll, data = LGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = LGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = LGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = LGSH)

kruskal.test(Proximity_to_FM ~ Year, data = LGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = LGSH)
dunnTest(Proximity_to_FM ~ Year, data = LGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = LGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = LGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = LGSH, method="bonferroni")

# individual species
# cytherea
kruskal.test(Proximity_to_FM ~ Year, data = cythereaLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = cythereaLGSH)
dunnTest(Proximity_to_FM ~ Year, data = cythereaLGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = cythereaLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = cythereaLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = cythereaLGSH, method="bonferroni")


# digitifera
kruskal.test(Proximity_to_FM ~ Atoll, data = digitiferaLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = digitiferaLGSH)


# humilis
kruskal.test(Proximity_to_FM ~ Atoll, data = humilisLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = humilisLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = humilisLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = humilisLGSH)

kruskal.test(Proximity_to_FM ~ Year, data = humilisLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = humilisLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = humilisLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = humilisLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = humilisLGSH, method="bonferroni")


# millepora
kruskal.test(Proximity_to_FM ~ Atoll, data = milleporaLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = milleporaLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = milleporaLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = milleporaLGSH)

kruskal.test(Proximity_to_FM ~ Year, data = milleporaLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = milleporaLGSH)
dunnTest(Proximity_to_FM ~ Year, data = milleporaLGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = milleporaLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = milleporaLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = milleporaLGSH, method="bonferroni")


# muricata
kruskal.test(Proximity_to_FM ~ Atoll, data = muricataLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = muricataLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = muricataLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = muricataLGSH)


# plantaginea
kruskal.test(Proximity_to_FM ~ Atoll, data = plantaLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = plantaLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = plantaLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = plantaLGSH)

kruskal.test(Proximity_to_FM ~ Year, data = plantaLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = plantaLGSH)
dunnTest(Proximity_to_FM ~ Year, data = plantaLGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = plantaLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = plantaLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = plantaLGSH, method="bonferroni")


# secale
kruskal.test(Proximity_to_FM ~ Atoll, data = secaleLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = secaleLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = secaleLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = secaleLGSH)

kruskal.test(Proximity_to_FM ~ Year, data = secaleLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = secaleLGSH)
dunnTest(Proximity_to_FM ~ Year, data = secaleLGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = secaleLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = secaleLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = secaleLGSH, method="bonferroni")


# tenuis
kruskal.test(Proximity_to_FM ~ Atoll, data = tenuisLGSH)
kruskal_effsize(Proximity_to_FM ~ Atoll, data = tenuisLGSH)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = tenuisLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Atoll, data = tenuisLGSH)

kruskal.test(Proximity_to_FM ~ Year, data = tenuisLGSH)
kruskal_effsize(Proximity_to_FM ~ Year, data = tenuisLGSH)
dunnTest(Proximity_to_FM ~ Year, data = tenuisLGSH, method="bonferroni")

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Year, data = tenuisLGSH)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Year, data = tenuisLGSH)
dunnTest(Spawning_Time_After_Sunset_Minutes ~ Year, data = tenuisLGSH, method="bonferroni")



#### wild v frame ####
# humilis, digi, and milli only, BAA only
# subset by species & Baa : due to sample size
humBaa <- subset(humilisLGSH, Atoll=="Baa")
digiBaa <- subset(digitiferaLGSH, Atoll=="Baa")
milliBaa <- subset(milleporaLGSH, Atoll=="Baa")
# all only 2 reef types --> frame & wild

kruskal.test(Proximity_to_FM ~ Reef_type, data = humBaa)
kruskal_effsize(Proximity_to_FM ~ Reef_type, data = humBaa)

kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Reef_type, data = humBaa)
kruskal_effsize(Spawning_Time_After_Sunset_Minutes ~ Reef_type, data = humBaa)


kruskal.test(Proximity_to_FM ~ Reef_type, data = digiBaa)
kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Reef_type, data = digiBaa)



kruskal.test(Proximity_to_FM ~ Reef_type, data = milliBaa)
kruskal.test(Spawning_Time_After_Sunset_Minutes ~ Reef_type, data = milliBaa)

# all not significant







#### Environmental cues ####

# Spearman's rank-order correlation
correlations <- LGSH %>%
  select(Proximity_to_FM, WindSpeed_DailyAV, Total_Daily_Precipitation_mm, Tide_Depth_m) %>%
  cor(method = "spearman", use = "complete.obs")
print(correlations)

#  Kendall's Tau-b correlation test 
# Load the necessary library (if not already loaded)
library(Kendall)
EnvLGSHNM <- subset(LGSH, Atoll=="NorthMale")
EnvLGSHBaa <- subset(LGSH, Atoll=="Baa")

# Correlation tests: LGSH
cor.test(x = LGSH$WindSpeed_DailyAV, y = LGSH$Proximity_to_FM, method = "kendall")
cor.test(x = LGSH$Total_Daily_Precipitation_mm, y = LGSH$Proximity_to_FM, method = "kendall")
cor.test(x = LGSH$Tide_Depth_m, y = LGSH$Proximity_to_FM, method = "kendall")
cor.test(x = LGSH$AverageDaily_SST, y = LGSH$Proximity_to_FM, method = "kendall")

# Correlation tests: Baa
cor.test(x = EnvLGSHBaa$WindSpeed_DailyAV, y = EnvLGSHBaa$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHBaa$Total_Daily_Precipitation_mm, y = EnvLGSHBaa$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHBaa$Tide_Depth_m, y = EnvLGSHBaa$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHBaa$AverageDaily_SST, y = EnvLGSHBaa$Proximity_to_FM, method = "kendall")

# Correlation tests: NM
cor.test(x = EnvLGSHNM$WindSpeed_DailyAV, y = EnvLGSHNM$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHNM$Total_Daily_Precipitation_mm, y = EnvLGSHNM$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHNM$Tide_Depth_m, y = EnvLGSHNM$Proximity_to_FM, method = "kendall")
cor.test(x = EnvLGSHNM$AverageDaily_SST, y = EnvLGSHNM$Proximity_to_FM, method = "kendall")

# data frame all
ppt <- c(0.2763633, 0.1154983, 0.4185969)
wind <- c(0.4811039, 0.2459261, 0.4377158)
tide <- c(0.6257439, 0.7243457, 0.4703455)
atolls <- c("both", "Baa", "NM")
kendall <- data.frame(atolls,ppt,wind,tide)

# data frame atolls only
ppt2 <- c(0.1154983, 0.4185969)
wind2 <- c(0.2459261, 0.4377158)
tide2 <- c(0.7243457, 0.4703455)
atolls2 <- c("Baa", "NM")
kendall2 <- data.frame(atolls2,ppt2,wind2,tide2)

# low values but p < 0.05: some degree of association between the variables, albeit not a strong one