#-------------------download and install these libraries
library(summarytools)
library(tidyverse)
library(readr)
library(readxl)
library(ggstatsplot)
library(lubridate)
library(RColorBrewer)
library(car)
library(dlookr)
library(patchwork)
library(statsExpressions)
library(egg)

#-------------------Set working directory and load datafile
setwd("D:/Proiect_TE/activitate_nocturna/Data&Analyses") # change to your own working directory
nocturnal_obs <-
  read_csv("Supplemental Data S2.csv")

#-------------------select only nocturnal observations from inat
nocturnal_inat <- filter(nocturnal_obs, source == "inat")

#-------------------add a year column and a month one
nocturnal_inat$year <- year(ymd(nocturnal_inat$date))
nocturnal_inat$month <- month(ymd(nocturnal_inat$date), label = T)

#-------------------import all inat data retrieved by 23rd Feb 2023
natrix_inat_t <- read_xlsx("Supplemental Data S1.xlsx")

#-------------------obtain date for every record
natrix_inat_t2 <- as_tibble(natrix_inat_t) %>%
  dplyr::select(observed_on) %>%
  rename("date" = "observed_on")

#-------------------make a vector for all data
type1 <- "all"

#-------------------add the vector for every observation
natrix_inat_t2$year <- year(ymd(natrix_inat_t2$date))
natrix_inat_t3 <- data.frame(type = rep(type1, each = 9745))
natrix_inat_t2$type <- natrix_inat_t3$type

#-------------------prepare data for inat bar plot
nocturnal_inat2 <- nocturnal_inat %>%
  dplyr::select(date, year, test_column) %>%
  mutate(
    type = case_when(
      test_column == "crepuscular" ~ "nocturnal",
      test_column == "nocturnal" ~ "nocturnal"
    )
  ) %>%
  dplyr::select(date, year, type)

#-------------------put nocturnal observations together with all inat data
natrix_d_combined <- rbind(natrix_inat_t2, nocturnal_inat2)
natrix_d_combined2 <- natrix_d_combined %>%
  dplyr::select(year, type) %>%
  group_by(year, type) %>%
  summarise(value = n()) %>%
  na.omit()

#-------------------make figure 2A
inat_by_year <- subset(natrix_d_combined2,
                       year %in% c(2014, 2015,
                                   2016, 2017, 2018, 2019, 2020, 2021, 2022)) %>%
  ggplot(aes(x = year, y = value, fill = type)) +
  geom_bar(stat = "identity",
           position = position_stack()) +
  geom_text(aes(label = value),
            fontface = "bold",
            vjust = -1) +
  labs(x = NULL, y = "\n Number of observations") +
  scale_fill_brewer(palette = "Accent") +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    legend.text = element_text(face = "bold",
                               margin = margin(r = 30, unit = "pt")),
    legend.position = "bottom",
    axis.title.x = element_text(face = "bold",
                                colour = "black"),
    axis.title.y = element_text(face = "bold",
                                colour = "black"),
    axis.text.x = element_text(face = "bold",
                               colour = "black"),
    axis.text.y = element_text(face = "bold",
                               colour = "black")
  )

#-------------------prepare data for figure 2b
nocturnal_obs$month <- month(ymd(nocturnal_obs$date), label = T)
nocturnal_obs_p <- nocturnal_obs %>%
  dplyr::select(`Body pattern`,
                test_column,
                ind_nr,
                `Age class`,
                `Medium (Aquatic or Terrestrial)`,
                month)

#-------------------make figure 2b
n_all_month <- nocturnal_obs_p %>%
  dplyr::select(month, test_column, ind_nr) %>%
  group_by(month, test_column) %>%
  summarise(ind_nr = sum(ind_nr)) %>%
  ggplot(aes(month, ind_nr, fill = test_column)) +
  geom_bar(
    stat = "identity",
    position = position_dodge2(preserve = "single"),
    alpha = .75
  ) +
  geom_text(
    aes(label = ind_nr),
    fontface = "bold",
    vjust = -1,
    position = position_dodge(.9)
  ) +
  labs(x = NULL, y = "\n Number of observed snakes") +
  scale_fill_brewer(palette = "Accent") +
  theme_classic() +
  theme(
    legend.title = element_blank(),
    legend.text = element_text(face = "bold",
                               margin = margin(r = 30, unit = "pt")),
    legend.position = "bottom",
    axis.title.x = element_text(face = "bold",
                                colour = "black"),
    axis.title.y = element_text(face = "bold",
                                colour = "black"),
    axis.text.x = element_text(face = "bold",
                               colour = "black"),
    axis.text.y = element_text(face = "bold",
                               colour = "black")
  )

#-------------------plot the 2 figs together
par(mar = c(4, 4, 0.1, 0.1))
paired_p <-
  ggarrange(
    inat_by_year,
    n_all_month,
    ncol = 2,
    labels = c("(A)", "(B)"),
    label.args = list(gp = grid::gpar(fontsize = 20))
  )
#save them as pdf
pdf(
  file = "Figure2.pdf",
  width = 11,
  height = 8.5,
  onefile = T
)
paired_p
dev.off()

#-------------------read data for violin plots and chi-squared
nocturnal_obs2 <-
  read_csv("Supplemental Data S2.csv")
nocturnal_obs2$moon_act <-
  ifelse(nocturnal_obs2$moon_up < 1, "no", "yes")

#-------------------assign to a new dataset to work with
nocturnal_morph <- nocturnal_obs2

#-------------------check normality of data
nocturnal_morph2 <- nocturnal_morph %>%
  dplyr::select(`Body pattern`, mean_temp_2m) %>%
  group_by(`Body pattern`) %>%
  normality(mean_temp_2m)
view(nocturnal_morph2)

#-------------------check for equal variance
leveneTest(mean_temp_2m ~ `Body pattern`, nocturnal_morph)

#-------------------set the y axis title
ylab_t <-
  expression(Mean ~ temperature ~ (T[24 * h]) ~ "in" ~ C * degree)

#-------------------realize the ANOVA with all groups
plot_morph <- ggbetweenstats(
  nocturnal_morph,
  `Body pattern`,
  mean_temp_2m,
  type = "p",
  plot.type = "box",
  var.equal = T,
  results.subtitle = T,
  xlab = NULL,
  ylab = ylab_t,
  bf.message = F,
  ggsignif.args = list(tip_length = 0.01,
                       na.rm = TRUE),
  ggtheme = ggplot2::theme_classic(),
  ggplot.component = list(
    theme(
      legend.position = "bottom",
      axis.title.x = element_blank(),
      axis.text = element_text(colour = "black")
    )
  ),
  centrality.label.args = list(alpha = 0, size = 5)
)

#-------------------check the statistic of the anova
extract_stats(plot_morph)

#-------------------filter only observation from sympatric region
nocturnal_morph_s <- nocturnal_obs2 %>% 
  filter(between(lat, 43, 47)) %>% 
  mutate(`Body pattern` = ifelse(`Body pattern` == "unstriped",
                                 "common morph", `Body pattern`)) %>% 
  drop_na(`Body pattern`)

#-------------------check normality of data
nocturnal_morph_s %>%
  dplyr::select(`Body pattern`, mean_temp_2m) %>%
  group_by(`Body pattern`) %>%
  normality(mean_temp_2m)

#-------------------check for equal variance
leveneTest(mean_temp_2m ~ `Body pattern`, nocturnal_morph_s)

#-------------------realize the K-W test with all groups
plot_morph_s <- ggbetweenstats(
  nocturnal_morph_s,
  `Body pattern`,
  mean_temp_2m,
  type = "n",
  plot.type = "box",
  var.equal = T,
  results.subtitle = T,
  xlab = NULL,
  ylab = ylab_t,
  bf.message = F,
  ggsignif.args = list(tip_length = 0.01,
                       na.rm = TRUE),
  violin.args = list(width = 0, alpha = 0, na.rm = TRUE),
  ggtheme = ggplot2::theme_classic(),
  ggplot.component = list(
    theme(
      legend.position = "bottom",
      axis.title.x = element_blank(),
      axis.text = element_text(colour = "black")
    )
  ),
  centrality.label.args = list(alpha = 0, size = 5)
)

#-------------------check the statistic of the K-W test
extract_stats(plot_morph)

#-------------------save it as pdf
pdf(
  file = "Figure4.pdf",
  width = 11,
  height = 8.5,
  onefile = T
)
plot_morph_s
dev.off()

#-------------------prepare data for mann-whitney test of T24h and type of behavior
nocturnal_obs3 <- filter(nocturnal_obs2, !is.na(test_column))

#-------------------check normality of data
nocturnal_obs4 <- nocturnal_obs3 %>%
  dplyr::select(test_column, mean_temp_2m) %>%
  group_by(test_column) %>%
  normality(mean_temp_2m)
view(nocturnal_obs4)

#-------------------check for equal variance
leveneTest(mean_temp_2m ~ test_column, nocturnal_obs3)

#-------------------perform mann-whitney
mann_t24 <- ggbetweenstats(nocturnal_obs3,
                           test_column,
                           mean_temp_2m,
                           type = "n")

testing <-
  wilcox.test(mean_temp_2m ~ test_column,
              nocturnal_obs3,
              exact = F,
              conf.int = T)
testing

#-------------------obtain a z-score
za <- qnorm(testing$p.value / 2)
za

#-------------------prepare data for mann-whitney test of moonlight fraction
#-------------------and type of behavior
nocturnal_obs5 <- nocturnal_obs3 %>%
  dplyr::filter(moon_act == "yes")

#-------------------check normality of data
nocturnal_obs5 %>%
  dplyr::select(test_column, moon_fraction) %>%
  group_by(test_column) %>%
  normality(moon_fraction)

#-------------------check for equal variance
leveneTest(moon_fraction ~ test_column, nocturnal_obs5)

#-------------------perform mann-whitney
mann_moon <- ggbetweenstats(nocturnal_obs5,
                            test_column,
                            moon_fraction,
                            type = "n")

another_test <-
  wilcox.test(moon_fraction ~ test_column,
              nocturnal_obs5,
              exact = F,
              conf.int = T)
another_test

#-------------------obtain a z-score
za2 <- qnorm(another_test$p.value / 2)
za2

#-------------------prepare data for chi-square test of moon presence and Age class
nocturnal_moon <-
  filter(nocturnal_morph, `Age class` != "U", !is.na(`Body pattern`))

#-------------------perform chi-square test
chi_age <- ggbarstats(
  data = nocturnal_moon,
  x = moon_act,
  k = 3,
  y = `Age class`,
  label = "BOTH"
)
chi_age
#-------------------extract stats of chi-square
chi_age_stat <- extract_stats(chi_age)
chi_age_stat

#-------------------perform chi-square test of moon presence and behavior
chi_behavior <- ggbarstats(data = nocturnal_morph,
                           x = moon_act,
                           y = test_column,
                           label = "BOTH")
chi_behavior

#-------------------extract stats of chi-square
chi_behavior_stat <- extract_stats(chi_behavior)
chi_behavior_stat