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

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

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

## Plot and Statistical Analysis

# Plot for Fig. 4
data_PulseRate_grp <- data_PulseRate %>% group_by(Species,Pulse_rate)
summary_PulseRate <- data_PulseRate_grp %>% summarize(Mean = mean(Normalized_responses), N =n(), CI = 1.96*(sd(Normalized_responses)/sqrt(N)))
summary_PulseRate$Species <- factor(summary_PulseRate$Species, levels = c("H.versicolor","H.chrysoscelis"))

ggplot(summary_PulseRate, aes(x = log10(Pulse_rate), y =Mean, shape = Species, fill= Species)) +
  geom_point(aes(y =Mean),size = 4) + geom_line()+
  geom_errorbar(width=.02, aes(ymin=Mean-CI, ymax=Mean+CI)) +
  scale_shape_manual(values = c(1,19)) + scale_fill_manual(values=c("#FFFFFF","#CCCCCC")) +
  theme_classic() + ylab("Normalized responses (%)") +
  scale_x_continuous("Pulse rate",breaks = c(log10(5), log10(8), log10(10), log10(20),
                                             log10(30), log10(40), log10(50), log10(60), log10(80),
                                             log10(100), log10(125)), 
                     labels = paste0(c(5,8,10,20,30,40,50,60,80,100,125))) +
  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. 4

# For H. chrysoscelis subjecs
data_PulseRate_Hch <- data_PulseRate %>% filter(Species == "H.chrysoscelis")
data_PulseRate_Hch$Pulse_rate <- factor(data_PulseRate_Hch$Pulse_rate, levels=c("50","20","30","40","60","80","125"))

lmer.hch_pr1 <- lmer(Normalized_responses ~ 1+ Pulse_rate + (1|Frog_ID), data=data_PulseRate_Hch, REML = FALSE)
lmer.hch_pr2 <- lmer(Normalized_responses ~ 1 + (1|Frog_ID), data=data_PulseRate_Hch, REML = FALSE)
anova(lmer.hch_pr1,lmer.hch_pr2)
summary(lmer.hch_pr1)

# For H.versicolor subjects
data_PulseRate_Hv <- data_PulseRate %>% filter(Species == "H.versicolor")
data_PulseRate_Hv$Pulse_rate <- factor(data_PulseRate_Hv$Pulse_rate, levels=c("20","5","8","10","30","50","100"))

lmer.hv_pr1 <- lmer(Normalized_responses ~ 1+ Pulse_rate + (1|Frog_ID), data=data_PulseRate_Hv, REML = FALSE)
lmer.hv_pr2 <- lmer(Normalized_responses ~ 1 + (1|Frog_ID), data=data_PulseRate_Hv, REML = FALSE)
anova(lmer.hv_pr1,lmer.hv_pr2)
summary(lmer.hv_pr1)