# Load required libraries
pacman::p_load(pacman, tidyverse, tidyr, dplyr, readr, grid, gridExtra, ggplot2, RColorBrewer)

# Load the CSV file 
data <- read.csv("Ferraro_PhanChorus_Sound.csv")

# Convert Trail and Treatment to factors for proper grouping
data$Trail <- as.factor(data$Trail)
data$Treatment <- as.factor(data$Treatment)

# Subtitle
Subtitle <- "Ordinal Date: 52 days - Jul 15th to Sept 4th; Treatment: “T” = phantom chorus was on. “C” = phantom chorus playback was off."
Subtitle1 <- "Treatment: “T” = phantom chorus was on. “C” = phantom chorus playback was off."

# Viz 1
# Aggregate data to calculate mean dbAL50 by Trail, Treatment, and Hour
daily_pattern <- data %>%
  group_by(Trail, Treatment, Hour) %>%
  summarize(mean_dbAL50 = mean(dbAL50, na.rm = TRUE))
## `summarise()` has grouped output by 'Trail', 'Treatment'. You can override
## using the `.groups` argument.
# Create the plot
ggplot(daily_pattern, aes(x = Hour, y = mean_dbAL50, color = Treatment)) +
  geom_line(size = 1) +  # Thicker lines for greater visibility
  facet_wrap(~ Trail) +  # Separate panels for each trail
  labs(title = "Average Daily Sound Level Pattern by Treatment and Trail",
       subtitle = Subtitle1,
       x = "Hour of Day",
       y = "Mean dbAL50 (dB)") +
  theme_bw() +  # Clean, white background theme
  scale_color_brewer(palette = "Set1")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Viz 2
# Aggregate data to calculate mean dbAL50 by Trail, Treatment, and OrdinalDate
trend_over_days <- data %>%
  group_by(Trail, Treatment, OrdinalDate) %>%
  summarize(mean_dbAL50 = mean(dbAL50, na.rm = TRUE))
## `summarise()` has grouped output by 'Trail', 'Treatment'. You can override
## using the `.groups` argument.
# Create the plot
ggplot(trend_over_days, aes(x = OrdinalDate, y = mean_dbAL50, color = Treatment)) +
  geom_line(size = 1) +
  facet_wrap(~ Trail) +
  labs(title = "Average Sound Level Over Days by Treatment and Trail",
       subtitle = Subtitle,
       x = "Ordinal Date",
       y = "Mean dbAL50 (dB)") +
  theme_bw() +
  scale_color_brewer(palette = "Set1")

# Viz 3
# Create the plot with violin and boxplot layers
ggplot(data, aes(x = Treatment, y = dbAL50, fill = Treatment)) +
  geom_violin(alpha = 0.7) +  # Semi-transparent violins for distribution shape
  geom_boxplot(width = 0.1, fill = "white") +  # Narrow boxplots inside violins
  facet_wrap(~ Trail) +
  labs(title = "Distribution of Sound Levels by Treatment and Trail",
       subtitle = Subtitle1,
       x = "Treatment",
       y = "dbAL50 (dB)") +
  theme_bw() +
  scale_fill_brewer(palette = "Set1")