## This script analyzes the normalized responses of H.chrysoscelis subjects to conspecific and heterospecific calls.  
# By running this script, you will be able to make figures used in Figure3 and perform statistical analysis associated with this figure.

# Read packages
library("dplyr")
library("ggplot2")
library("lme4")
library("lmerTest")
library("binom")

## Read data 
data_SpDisc <- read.csv("Figure3_SpeciesDiscrimination.csv",header = T, check.names=F) # full data

## Plot and Statistical analysis 

# Plot for Fig 3a (Breeding season)
data_SpDisc_breed <- data_SpDisc %>% filter(Season == "Breeding")
data_SpDisc_breed_grp <- data_SpDisc_breed %>% group_by(Stimulus, Amplitude)
summary_SpDisc_breed <- data_SpDisc_breed_grp %>% summarise(Mean = mean(Normalized_responses), N =n(), CI = 1.96*(sd(Normalized_responses)/sqrt(N)))
summary_SpDisc_breed$Amplitude <- factor(summary_SpDisc_breed$Amplitude, c(65,75,85))
summary_SpDisc_breed$Stimulus <- factor(summary_SpDisc_breed$Stimulus, c("H.versicolor","H.chrysoscelis"))
ggplot(summary_SpDisc_breed, aes(x=Amplitude,y=Mean, fill=Stimulus)) +
  geom_bar(position=position_dodge(.9),width=.6, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.2, aes(ymin=Mean-CI, ymax=Mean+CI)) +
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC")) + theme_classic() + geom_hline(yintercept=0, linetype=5) +
  ylim(-20,150) + ylab("Normalized responses(%)") + xlab("Signal level (dB SPL)") +
  theme(axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14),axis.text.y=element_text(size=14),legend.title=element_blank(),legend.position="top",panel.spacing.x=unit(1, "lines"))

# Analysis for Fig. 3a
data_SpDisc_breed$Amplitude <- factor(data_SpDisc_breed$Amplitude, levels = c(85,65,75))
data_SpDisc_breed$Stimulus <- factor(data_SpDisc_breed$Stimulus)
lmer.breed1 <- lmer(Normalized_responses ~ Amplitude + Stimulus + Amplitude*Stimulus + (1|Frog_ID), data=data_SpDisc_breed, REML = FALSE)
lmer.breed2 <- lmer(Normalized_responses ~ Amplitude + Stimulus + (1|Frog_ID), data=data_SpDisc_breed, REML = FALSE)
anova(lmer.breed1,lmer.breed2)
summary(lmer.breed1)

# For Fig 3b (Non-breeding season)
data_SpDisc_nobreed <- data_SpDisc %>% filter(Season == "Non-breeding")
data_SpDisc_nobreed_grp <- data_SpDisc_nobreed %>% group_by(Stimulus, Amplitude)
summary_SpDisc_nobreed <- data_SpDisc_nobreed_grp %>% summarise(Mean = mean(Normalized_responses), N =n(), CI = 1.96*(sd(Normalized_responses)/sqrt(N)))
summary_SpDisc_nobreed$Amplitude <- factor(summary_SpDisc_nobreed$Amplitude, c(85))
summary_SpDisc_nobreed$Stimulus <- factor(summary_SpDisc_nobreed$Stimulus, c("H.versicolor","H.chrysoscelis"))
ggplot(summary_SpDisc_nobreed, aes(x=Amplitude,y=Mean, fill=Stimulus)) +
  geom_bar(position=position_dodge(.9),width=.6, colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.2, aes(ymin=Mean-CI, ymax=Mean+CI)) +
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC")) + theme_classic() + geom_hline(yintercept=0, linetype=5) +
  ylim(-20,150) + ylab("Normalized responses(%)") + xlab("Signal level (dB SPL)") +
  theme(axis.title.x=element_text(size=16), axis.title.y=element_text(size=16), axis.text.x=element_text(size=14),axis.text.y=element_text(size=14),legend.title=element_blank(),legend.position="top",panel.spacing.x=unit(1, "lines"))

# Analysis for Fig. 3b
data_SpDisc_nobreed$Amplitude <- factor(data_SpDisc_nobreed$Amplitude)
data_SpDisc_nobreed$Stimulus <- factor(data_SpDisc_nobreed$Stimulus)
lmer.nobreed <- lmer(Normalized_responses ~ Stimulus + (1|Frog_ID), data=data_SpDisc_nobreed, REML = FALSE)
summary(lmer.nobreed)

# Analysis on combined data (breeding + non-breeding season)
data_SpDisc$Amplitude <- factor(data_SpDisc$Amplitude)
data_SpDisc$Stimulus <- factor(data_SpDisc$Stimulus)

data_SpDisc_85 <- data_SpDisc %>% filter(Amplitude == 85) # Filtering 85dB to compare data from the two seasons

data_SpDisc_85_cons <- data_SpDisc_85 %>% filter(Stimulus == "H.chrysoscelis") # Comparing only conspecific sounds
lmer.1 <- lm(Normalized_responses ~ 1+ Season, data=data_SpDisc_85_cons)
summary(lmer.1)