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

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

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

## Plot and Statistical Analysis

# Plot and psychometric curve in Fig. 6
data_CallAmp_grp <- data_CallAmp %>% group_by(Amplitude)
summary_CallAmp_grp <- data_CallAmp_grp %>% summarize(Mean =mean(Normalized_responses), N = n(), CI = 1.96*(sd(Normalized_responses)/sqrt(N)))

Amplitude <- summary_CallAmp_grp$Amplitude
Mean <- summary_CallAmp_grp$Mean
df_ref <- data.frame(Amplitude,Mean) # Making a data frame of amplitude and normalized responses

p <- nls(Mean ~ SSlogis(Amplitude, Asym, xmid, scal), data = df_ref)
pframe <- data.frame(Amplitude=10^seq(log10(0+1),log10(85), length.out=100))
pp <- data.frame(pframe, Mean = predict(p,pframe), wts =FALSE) # psychometric data

ggplot(summary_CallAmp_grp, aes(x = Amplitude, y = Mean)) +
  geom_point(aes(y =Mean),size = 4, shape = 16) +
  geom_errorbar(width=.7, aes(ymin=Mean-CI, ymax=Mean+CI)) +
  theme_classic() + ylab("Normalized responses (%)") +
  scale_x_continuous(breaks = c(0,25,31,37,43,49,55,61,67,73,79,85), 
                     labels = paste0(c(0,25,31,37,43,49,55,61,67,73,79,85))) +
  geom_line(data =pp, aes(linetype=wts), size=1) + scale_y_continuous(breaks=seq(-25,150,25)) +
  geom_hline(yintercept=0, linetype=5) +
  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="none",panel.spacing.x=unit(1, "lines"))

# Statistical analysis for Fig. 6
data_CallAmp$Amplitude <- factor(data_CallAmp$Amplitude, levels = c(0,25,31,37,43,49,55,61,67,73,79,85))

lmer.1 <- lmer(Normalized_responses ~ 1 + Amplitude + (1|Frog_ID), data_CallAmp, REML = FALSE)
lmer.2 <- lmer(Normalized_responses ~ 1  + (1|Frog_ID), data_CallAmp, REML = FALSE)
anova(lmer.1,lmer.2)
summary(lmer.1)