## This script analyzes choice data collected for Experiment 2.  
# By running this script, you will be able to make figures used in Figure 4and5 and perform statistical analysis associated with these figures

# Read packages
library("dplyr")
library("ggplot2")
library("binom")
library("geepack")
library("emmeans")

## Read data for Exp 2 
data_nochoice_1.25kHz <- read.csv("Figure4&5_SignalRecognition_1.25kHz_IM.csv",header = T, check.names=F) # full data
data_nochoice_1.25kHz$ID = factor(data_nochoice_1.25kHz$ID)
data_nochoice_2.5kHz <- read.csv("Figure4&5_SignalRecognition_2.5kHz_IM.csv",header = T, check.names=F) # full data
data_nochoice_2.5kHz$ID = factor(data_nochoice_2.5kHz$ID)

# Analyze data for response rate when target signal is 1.25 kHz

  # Pooled across the three TMRs
  data_nochoice_1.25kHz_pooled <- data_nochoice_1.25kHz %>% group_by(Condition)
  summary_nochoice_1.25kHz_pooled <- data_nochoice_1.25kHz_pooled %>% summarise(Response = sum(Response), N =n())
  binomial_nochoice_1.25kHz_pooled <- binom.exact(summary_nochoice_1.25kHz_pooled$Response, summary_nochoice_1.25kHz_pooled$N) %>%
    mutate(Condition = summary_nochoice_1.25kHz_pooled$Condition)
  binomial_nochoice_1.25kHz_pooled$Condition = factor(binomial_nochoice_1.25kHz_pooled$Condition,c("Signal alone","IM Test","EM Control"))
  binomial_nochoice_1.25kHz_pooled$mean <- binomial_nochoice_1.25kHz_pooled$mean *100 # Convert proportions to percentage responding
  binomial_nochoice_1.25kHz_pooled$lower <- binomial_nochoice_1.25kHz_pooled$lower *100 # Convert proportions to percentage responding
  binomial_nochoice_1.25kHz_pooled$upper <- binomial_nochoice_1.25kHz_pooled$upper *100 # Convert proportions to percentage responding
  
  # For each TMR
  data_nochoice_1.25kHz_TMR <- data_nochoice_1.25kHz %>% group_by(TMR,Condition)
  summary_nochoice_1.25kHz_TMR <- data_nochoice_1.25kHz_TMR %>% summarise(Response = sum(Response), N =n())
  binomial_nochoice_1.25kHz_TMR <- binom.exact(summary_nochoice_1.25kHz_TMR$Response, summary_nochoice_1.25kHz_TMR$N) %>%
    mutate(TMR = summary_nochoice_1.25kHz_TMR$TMR) %>% mutate(Condition = summary_nochoice_1.25kHz_TMR$Condition)
  binomial_nochoice_1.25kHz_TMR$TMR = factor(binomial_nochoice_1.25kHz_TMR$TMR,c("0dB","-6dB","-12dB"))
  binomial_nochoice_1.25kHz_TMR$Condition = factor(binomial_nochoice_1.25kHz_TMR$Condition,c("Signal alone","IM Test","EM Control"))
  binomial_nochoice_1.25kHz_TMR$mean <- binomial_nochoice_1.25kHz_TMR$mean *100 # Convert proportions to percentage responding
  binomial_nochoice_1.25kHz_TMR$lower <- binomial_nochoice_1.25kHz_TMR$lower *100 # Convert proportions to percentage responding
  binomial_nochoice_1.25kHz_TMR$upper <- binomial_nochoice_1.25kHz_TMR$upper *100 # Convert proportions to percentage responding
  
  
# Analyze threshold and latency data
data_nochoice_1.25kHz_response <- data_nochoice_1.25kHz %>% filter(Response != 0)
data_nochoice_1.25kHz_response <- data_nochoice_1.25kHz_response %>% group_by(Condition, Trial_type, TMR)
summary_nochoice_1.25kHz_threshold <- data_nochoice_1.25kHz_response %>% summarise(Median_threshold = median(Threshold), quantile_1 =quantile(Threshold,0.25), quantile_3 =quantile(Threshold,0.75), Mean_threshold = mean(Threshold), N =n())
summary_nochoice_1.25kHz_latency <- data_nochoice_1.25kHz_response %>% summarise(Median_latency = median(Latency), quantile_1 =quantile(Latency,0.25), quantile_3 =quantile(Latency,0.75), Mean_latency = mean(Latency), N =n())
  
# Analyze data for response rate when target signal is 2.5kHz

# Pooled across the three TMRs
data_nochoice_2.5kHz_pooled <- data_nochoice_2.5kHz %>% group_by(Condition)
summary_nochoice_2.5kHz_pooled <- data_nochoice_2.5kHz_pooled %>% summarise(Response = sum(Response), N =n())
binomial_nochoice_2.5kHz_pooled <- binom.exact(summary_nochoice_2.5kHz_pooled$Response, summary_nochoice_2.5kHz_pooled$N) %>%
  mutate(Condition = summary_nochoice_2.5kHz_pooled$Condition)
binomial_nochoice_2.5kHz_pooled$Condition = factor(binomial_nochoice_2.5kHz_pooled$Condition,c("Signal alone","IM Test","EM Control"))
binomial_nochoice_2.5kHz_pooled$mean <- binomial_nochoice_2.5kHz_pooled$mean *100 # Convert proportions to percentage responding
binomial_nochoice_2.5kHz_pooled$lower <- binomial_nochoice_2.5kHz_pooled$lower *100 # Convert proportions to percentage responding
binomial_nochoice_2.5kHz_pooled$upper <- binomial_nochoice_2.5kHz_pooled$upper *100 # Convert proportions to percentage responding

# For each TMR
data_nochoice_2.5kHz_TMR <- data_nochoice_2.5kHz %>% group_by(TMR,Condition)
summary_nochoice_2.5kHz_TMR <- data_nochoice_2.5kHz_TMR %>% summarise(Response = sum(Response), N =n())
binomial_nochoice_2.5kHz_TMR <- binom.exact(summary_nochoice_2.5kHz_TMR$Response, summary_nochoice_2.5kHz_TMR$N) %>%
  mutate(TMR = summary_nochoice_2.5kHz_TMR$TMR) %>% mutate(Condition = summary_nochoice_2.5kHz_TMR$Condition)
binomial_nochoice_2.5kHz_TMR$TMR = factor(binomial_nochoice_2.5kHz_TMR$TMR,c("0dB","-6dB","-12dB"))
binomial_nochoice_2.5kHz_TMR$Condition = factor(binomial_nochoice_2.5kHz_TMR$Condition,c("Signal alone","IM Test","EM Control"))
binomial_nochoice_2.5kHz_TMR$mean <- binomial_nochoice_2.5kHz_TMR$mean *100 # Convert proportions to percentage responding
binomial_nochoice_2.5kHz_TMR$lower <- binomial_nochoice_2.5kHz_TMR$lower *100 # Convert proportions to percentage responding
binomial_nochoice_2.5kHz_TMR$upper <- binomial_nochoice_2.5kHz_TMR$upper *100 # Convert proportions to percentage responding


# Analyze threshold and latency data
data_nochoice_2.5kHz_response <- data_nochoice_2.5kHz %>% filter(Response != 0)
data_nochoice_2.5kHz_response <- data_nochoice_2.5kHz_response %>% group_by(Condition, Trial_type, TMR)
summary_nochoice_2.5kHz_threshold <- data_nochoice_2.5kHz_response %>% summarise(Median_threshold = median(Threshold), quantile_1 =quantile(Threshold,0.25), quantile_3 =quantile(Threshold,0.75), Mean_threshold = mean(Threshold), N =n())
summary_nochoice_2.5kHz_latency <- data_nochoice_2.5kHz_response %>% summarise(Median_latency = median(Latency), quantile_1 =quantile(Latency,0.25), quantile_3 =quantile(Latency,0.75), Mean_latency = mean(Latency), N =n())

## Plots for Fig 3b depicting response rates (pooled and for each TMR when target signal is 1.25kHz)

  # For Fig. 3B (response rate pooled across TMRs)
  data_nochoice_1.25kHz_response$Condition <- factor(data_nochoice_1.25kHz_response$Condition,c("Signal alone","IM Test","EM Control"))
  data_nochoice_1.25kHz_response$Trial_type <- factor(data_nochoice_1.25kHz_response$Trial_type,c("IM","EM"))
  ggplot(binomial_nochoice_1.25kHz_pooled, aes(x = Condition, y = mean, fill=Condition)) + 
    geom_bar(stat="identity",position=position_dodge()) + 
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position=position_dodge(.9)) + 
    theme_classic() + ylab("Response rate (%)") + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9")) + 
    ylim (0,100) + theme(axis.title.x=element_text(size=20),axis.title.y=element_text(size=20), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    xlab("Pooled")
  
  # For Fig. 3B (response rate for each TMR)
  ggplot(binomial_nochoice_1.25kHz_TMR, aes(x =TMR, y = mean, fill=Condition)) + 
    geom_bar(stat="identity",position=position_dodge(.9)) + 
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position=position_dodge(.9)) + 
    theme_classic() + ylab("Response rate (%)") + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9")) + 
    ylim (0,100) + theme_classic()  + theme(axis.title.x=element_text(size=20),axis.title.y=element_text(size=20), axis.text.x=element_text(size=18),axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines"))
  
## Plots for Fig 3d depicting response rates (pooled and for each TMR when target signal is 2.5kHz)

  # For Fig. 3D (response rate pooled across TMRs)
  data_nochoice_2.5kHz_response$Condition <- factor(data_nochoice_2.5kHz_response$Condition,c("Signal alone","IM Test","EM Control"))
  data_nochoice_2.5kHz_response$Trial_type <- factor(data_nochoice_2.5kHz_response$Trial_type,c("IM","EM"))
  ggplot(binomial_nochoice_2.5kHz_pooled, aes(x = Condition, y = mean, fill=Condition)) + 
    geom_bar(stat="identity",position=position_dodge()) + 
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position=position_dodge(.9)) + 
    theme_classic() + ylab("Response rate (%)") + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9")) + 
    ylim (0,100) + theme(axis.title.x=element_text(size=20),axis.title.y=element_text(size=20), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    xlab("Pooled")
  
  # For Fig. 3D (response rate for each TMR)
  ggplot(binomial_nochoice_2.5kHz_TMR, aes(x =TMR, y = mean, fill=Condition)) + 
    geom_bar(stat="identity",position=position_dodge(.9)) + 
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.2, position=position_dodge(.9)) + 
    theme_classic() + ylab("Response rate (%)") + scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9")) + 
    ylim (0,100) + theme_classic()  + theme(axis.title.x=element_text(size=20),axis.title.y=element_text(size=20), axis.text.x=element_text(size=18),axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines"))
  
## Plots for Fig 4 
  
  # For Fig. 4B (depicting pulse number thresholds when target signal is 1.25kHz)
  source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")
  ggplot(data_nochoice_1.25kHz_response,aes(x=Trial_type, y=Threshold, fill=Condition)) + geom_flat_violin (position = position_nudge(x = .2, y = 0),alpha = .4) +
    geom_point(aes(y =Threshold, color = Condition), position=position_dodge(width = 0.3), size = 1, alpha = 0.4) + 
    geom_boxplot(width = 0.2,alpha=0.6, outlier.shape = NA, position=position_dodge(width = 0.3)) + scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
    scale_fill_manual(values =  c("#999999", "#E69F00", "#56B4E9")) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_text(size=20), axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    ylab("Pulse number threshold") + stat_summary(fun=mean, aes(group=Condition), geom="point", shape=23, size=2, col="black", fill="red", position=position_dodge(width = 0.3))+
    scale_y_continuous(limits = c(0,22), breaks = seq(0,20, by = 3)) + facet_grid(~TMR, scales = "fixed") + theme(aspect.ratio = 1)
  
  # For Fig. 4C (depicting response latencies when target signal is 1.25kHz)
  source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")
  ggplot(data_nochoice_1.25kHz_response,aes(x=Trial_type, y=Latency, fill=Condition)) + geom_flat_violin (position = position_nudge(x = .2, y = 0),alpha = .4) +
    geom_point(aes(y =Latency, color = Condition), position=position_dodge(width = 0.3), size = 1, alpha = 0.4) + 
    geom_boxplot(width = 0.2,alpha=0.6, outlier.shape = NA, position=position_dodge(width = 0.3)) + scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
    scale_fill_manual(values =  c("#999999", "#E69F00", "#56B4E9")) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_text(size=20), axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.text.y=element_text(size=16),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    ylab("Latency (ms)") + stat_summary(fun=mean, aes(group=Condition), geom="point", shape=23, size=2, col="black", fill="red", position=position_dodge(width = 0.3))+
    scale_y_continuous(limits = c(0,350), breaks = seq(0,300, by = 50)) + facet_grid(~TMR, scales = "fixed") + theme(aspect.ratio = 1)
  
  # For Fig. 4E (depicting pulse number thresholds when target signal is 1.25kHz)
  source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")
  ggplot(data_nochoice_2.5kHz_response,aes(x=Trial_type, y=Threshold, fill=Condition)) + geom_flat_violin (position = position_nudge(x = .2, y = 0),alpha = .4) +
    geom_point(aes(y =Threshold, color = Condition), position=position_dodge(width = 0.3), size = 1, alpha = 0.4) + 
    geom_boxplot(width = 0.2,alpha=0.6, outlier.shape = NA, position=position_dodge(width = 0.3)) + scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
    scale_fill_manual(values =  c("#999999", "#E69F00", "#56B4E9")) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_text(size=20), axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.text.y=element_text(size=18),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    ylab("Pulse number threshold") + stat_summary(fun=mean, aes(group=Condition), geom="point", shape=23, size=2, col="black", fill="red", position=position_dodge(width = 0.3))+
    scale_y_continuous(limits = c(0,22), breaks = seq(0,20, by = 3)) + facet_grid(~TMR, scales = "fixed") + theme(aspect.ratio = 1)
  
  # For Fig. 4F (depicting response latencies when target signal is 2.5kHz)
  source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")
  ggplot(data_nochoice_2.5kHz_response,aes(x=Trial_type, y=Latency, fill=Condition)) + geom_flat_violin (position = position_nudge(x = .2, y = 0),alpha = .4) +
    geom_point(aes(y =Latency, color = Condition), position=position_dodge(width = 0.3), size = 1, alpha = 0.4) + 
    geom_boxplot(width = 0.2,alpha=0.6, outlier.shape = NA, position=position_dodge(width = 0.3)) + scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
    scale_fill_manual(values =  c("#999999", "#E69F00", "#56B4E9")) + theme_classic() + theme(axis.title.x=element_blank(),axis.title.y=element_text(size=20), axis.text.x=element_blank(),axis.ticks.x=element_blank(), axis.text.y=element_text(size=16),legend.position="none",panel.spacing.x=unit(1, "lines")) +
    ylab("Latency (ms)") + stat_summary(fun=mean, aes(group=Condition), geom="point", shape=23, size=2, col="black", fill="red", position=position_dodge(width = 0.3))+
    scale_y_continuous(limits = c(0,350), breaks = seq(0,300, by = 50)) + facet_grid(~TMR, scales = "fixed") + theme(aspect.ratio = 1)
  
  
  ## Statistical analysis for Experiment 2 when target signals were 1.25 kHz
  
  # Using GEE to statistically analyze data for response rates
  data_nochoice_1.25kHz$Condition <- factor(data_nochoice_1.25kHz$Condition, c("Signal alone","IM Test","EM Control"))
  
    # Testing the effect of interaction between masking conditions and TMR
    gee.nochoice_1.25kHz_response.1  <- geeglm(Response ~  Condition + TMR + Condition:TMR, id = ID,
                                       data = data_nochoice_1.25kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
    gee.nochoice_1.25kHz_response.2  <- geeglm(Response ~  Condition + TMR, id = ID,
                                       data = data_nochoice_1.25kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
    anova(gee.nochoice_1.25kHz_response.1, gee.nochoice_1.25kHz_response.2)
    
    # Testing the effect of TMR
    gee.nochoice_1.25kHz_response.3  <- geeglm(Response ~  Condition, id = ID,
                                       data = data_nochoice_1.25kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
    anova(gee.nochoice_1.25kHz_response.2, gee.nochoice_1.25kHz_response.3)
    
    # Testing the effect of condition
    gee.nochoice_1.25kHz_response.4  <- geeglm(Response ~  1, id = ID,
                                               data = data_nochoice_1.25kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
    anova(gee.nochoice_1.25kHz_response.3, gee.nochoice_1.25kHz_response.4)
    
    # Testing the differences between the three conditions
    emmeans(gee.nochoice_1.25kHz_response.3, pairwise ~ Condition, adjust ="none")

  # Using Wilcoxon tests to statistically analyze data for thresholds and latencies
  
  data_nochoice_1.25kHz_IM_SA <- data_nochoice_1.25kHz %>% filter(Trial_type == "IM") %>% filter(Condition == "Signal alone")
  data_nochoice_1.25kHz_IM_IM <- data_nochoice_1.25kHz %>% filter(Trial_type == "IM") %>% filter(Condition == "IM Test")
  data_nochoice_1.25kHz_EM_SA <- data_nochoice_1.25kHz %>% filter(Trial_type == "EM") %>% filter(Condition == "Signal alone")
  data_nochoice_1.25kHz_EM_EM <- data_nochoice_1.25kHz %>% filter(Trial_type == "EM") %>% filter(Condition == "EM Control")
    
    # For TMR of -12dB
    data_nochoice_1.25kHz_IM_SA_m12dB <- data_nochoice_1.25kHz_IM_SA %>% filter(TMR == "-12dB")
    data_nochoice_1.25kHz_IM_IM_m12dB <- data_nochoice_1.25kHz_IM_IM %>% filter(TMR == "-12dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_1.25kHz_IM_SA_m12dB$Threshold, data_nochoice_1.25kHz_IM_IM_m12dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_IM_SA_m12dB$Latency, data_nochoice_1.25kHz_IM_IM_m12dB$Latency, paired = TRUE) 
    
    data_nochoice_1.25kHz_EM_SA_m12dB <- data_nochoice_1.25kHz_EM_SA %>% filter(TMR == "-12dB")
    data_nochoice_1.25kHz_EM_EM_m12dB <- data_nochoice_1.25kHz_EM_EM %>% filter(TMR == "-12dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_1.25kHz_EM_SA_m12dB$Threshold, data_nochoice_1.25kHz_EM_EM_m12dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_EM_SA_m12dB$Latency, data_nochoice_1.25kHz_EM_EM_m12dB$Latency, paired = TRUE) 
    
      # Difference between iM Test and EM Control
      wilcox.test(data_nochoice_1.25kHz_IM_IM_m12dB$Threshold, data_nochoice_1.25kHz_EM_EM_m12dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_1.25kHz_IM_IM_m12dB$Latency, data_nochoice_1.25kHz_EM_EM_m12dB$Latency, alternative = "two.sided")
      
    # For TMR of -6dB
      data_nochoice_1.25kHz_IM_SA_m6dB <- data_nochoice_1.25kHz_IM_SA %>% filter(TMR == "-6dB")
      data_nochoice_1.25kHz_IM_IM_m6dB <- data_nochoice_1.25kHz_IM_IM %>% filter(TMR == "-6dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_1.25kHz_IM_SA_m6dB$Threshold, data_nochoice_1.25kHz_IM_IM_m6dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_IM_SA_m6dB$Latency, data_nochoice_1.25kHz_IM_IM_m6dB$Latency, paired = TRUE) 
      
      data_nochoice_1.25kHz_EM_SA_m6dB <- data_nochoice_1.25kHz_EM_SA %>% filter(TMR == "-6dB")
      data_nochoice_1.25kHz_EM_EM_m6dB <- data_nochoice_1.25kHz_EM_EM %>% filter(TMR == "-6dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_1.25kHz_EM_SA_m6dB$Threshold, data_nochoice_1.25kHz_EM_EM_m6dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_EM_SA_m6dB$Latency, data_nochoice_1.25kHz_EM_EM_m6dB$Latency, paired = TRUE) 
      
      # Difference between IM Test and EM Control
      wilcox.test(data_nochoice_1.25kHz_IM_IM_m6dB$Threshold, data_nochoice_1.25kHz_EM_EM_m6dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_1.25kHz_IM_IM_m6dB$Latency, data_nochoice_1.25kHz_EM_EM_m6dB$Latency, alternative = "two.sided")
      
    # For TMR of 0dB
      data_nochoice_1.25kHz_IM_SA_0dB <- data_nochoice_1.25kHz_IM_SA %>% filter(TMR == "0dB")
      data_nochoice_1.25kHz_IM_IM_0dB <- data_nochoice_1.25kHz_IM_IM %>% filter(TMR == "0dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_1.25kHz_IM_SA_0dB$Threshold, data_nochoice_1.25kHz_IM_IM_0dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_IM_SA_0dB$Latency, data_nochoice_1.25kHz_IM_IM_0dB$Latency, paired = TRUE) 
      
      data_nochoice_1.25kHz_EM_SA_0dB <- data_nochoice_1.25kHz_EM_SA %>% filter(TMR == "0dB")
      data_nochoice_1.25kHz_EM_EM_0dB <- data_nochoice_1.25kHz_EM_EM %>% filter(TMR == "0dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_1.25kHz_EM_SA_0dB$Threshold, data_nochoice_1.25kHz_EM_EM_0dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_1.25kHz_EM_SA_0dB$Latency, data_nochoice_1.25kHz_EM_EM_0dB$Latency, paired = TRUE) 
      
      # Difference between IM Test and EM Control
      wilcox.test(data_nochoice_1.25kHz_IM_IM_0dB$Threshold, data_nochoice_1.25kHz_EM_EM_0dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_1.25kHz_IM_IM_0dB$Latency, data_nochoice_1.25kHz_EM_EM_0dB$Latency, alternative = "two.sided")
      
      
  ## Statistical analysis for Experiment 2 when target signals were 2.5 kHz
      # Using GEE to statistically analyze data for response rates
      data_nochoice_2.5kHz$Condition <- factor(data_nochoice_2.5kHz$Condition, c("Signal alone","IM Test","EM Control"))
      
      # Testing the effect of interaction between masking conditions and TMR
      gee.nochoice_2.5kHz_response.1  <- geeglm(Response ~  Condition + TMR + Condition:TMR, id = ID,
                                                data = data_nochoice_2.5kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
      gee.nochoice_2.5kHz_response.2  <- geeglm(Response ~  Condition + TMR, id = ID,
                                                data = data_nochoice_2.5kHz, family = binomial(link="logit"), corstr = "ex", std.err = "san.se")
      # Cannot fit glm model because of complete separation
      
      
      # Using Wilcoxon tests to statistically analyze data for thresholds and latencies
      
      data_nochoice_2.5kHz_IM_SA <- data_nochoice_2.5kHz %>% filter(Trial_type == "IM") %>% filter(Condition == "Signal alone")
      data_nochoice_2.5kHz_IM_IM <- data_nochoice_2.5kHz %>% filter(Trial_type == "IM") %>% filter(Condition == "IM Test")
      data_nochoice_2.5kHz_EM_SA <- data_nochoice_2.5kHz %>% filter(Trial_type == "EM") %>% filter(Condition == "Signal alone")
      data_nochoice_2.5kHz_EM_EM <- data_nochoice_2.5kHz %>% filter(Trial_type == "EM") %>% filter(Condition == "EM Control")
      
      # For TMR of -12dB
      data_nochoice_2.5kHz_IM_SA_m12dB <- data_nochoice_2.5kHz_IM_SA %>% filter(TMR == "-12dB")
      data_nochoice_2.5kHz_IM_IM_m12dB <- data_nochoice_2.5kHz_IM_IM %>% filter(TMR == "-12dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_2.5kHz_IM_SA_m12dB$Threshold, data_nochoice_2.5kHz_IM_IM_m12dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_IM_SA_m12dB$Latency, data_nochoice_2.5kHz_IM_IM_m12dB$Latency, paired = TRUE) 
      
      data_nochoice_2.5kHz_EM_SA_m12dB <- data_nochoice_2.5kHz_EM_SA %>% filter(TMR == "-12dB")
      data_nochoice_2.5kHz_EM_EM_m12dB <- data_nochoice_2.5kHz_EM_EM %>% filter(TMR == "-12dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_2.5kHz_EM_SA_m12dB$Threshold, data_nochoice_2.5kHz_EM_EM_m12dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_EM_SA_m12dB$Latency, data_nochoice_2.5kHz_EM_EM_m12dB$Latency, paired = TRUE) 
      
      # Difference between IM Test and EM Control
      wilcox.test(data_nochoice_2.5kHz_IM_IM_m12dB$Threshold, data_nochoice_2.5kHz_EM_EM_m12dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_2.5kHz_IM_IM_m12dB$Latency, data_nochoice_2.5kHz_EM_EM_m12dB$Latency, alternative = "two.sided")
      
      # For TMR of -6dB
      data_nochoice_2.5kHz_IM_SA_m6dB <- data_nochoice_2.5kHz_IM_SA %>% filter(TMR == "-6dB")
      data_nochoice_2.5kHz_IM_IM_m6dB <- data_nochoice_2.5kHz_IM_IM %>% filter(TMR == "-6dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_2.5kHz_IM_SA_m6dB$Threshold, data_nochoice_2.5kHz_IM_IM_m6dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_IM_SA_m6dB$Latency, data_nochoice_2.5kHz_IM_IM_m6dB$Latency, paired = TRUE) 
      
      data_nochoice_2.5kHz_EM_SA_m6dB <- data_nochoice_2.5kHz_EM_SA %>% filter(TMR == "-6dB")
      data_nochoice_2.5kHz_EM_EM_m6dB <- data_nochoice_2.5kHz_EM_EM %>% filter(TMR == "-6dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_2.5kHz_EM_SA_m6dB$Threshold, data_nochoice_2.5kHz_EM_EM_m6dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_EM_SA_m6dB$Latency, data_nochoice_2.5kHz_EM_EM_m6dB$Latency, paired = TRUE) 
      
      # Difference between IM Test and EM Control
      wilcox.test(data_nochoice_2.5kHz_IM_IM_m6dB$Threshold, data_nochoice_2.5kHz_EM_EM_m6dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_2.5kHz_IM_IM_m6dB$Latency, data_nochoice_2.5kHz_EM_EM_m6dB$Latency, alternative = "two.sided")
      
      # For TMR of 0dB
      data_nochoice_2.5kHz_IM_SA_0dB <- data_nochoice_2.5kHz_IM_SA %>% filter(TMR == "0dB")
      data_nochoice_2.5kHz_IM_IM_0dB <- data_nochoice_2.5kHz_IM_IM %>% filter(TMR == "0dB")
      # Difference between signal alone and IM Test
      wilcox.test(data_nochoice_2.5kHz_IM_SA_0dB$Threshold, data_nochoice_2.5kHz_IM_IM_0dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_IM_SA_0dB$Latency, data_nochoice_2.5kHz_IM_IM_0dB$Latency, paired = TRUE) 
      
      data_nochoice_2.5kHz_EM_SA_0dB <- data_nochoice_2.5kHz_EM_SA %>% filter(TMR == "0dB")
      data_nochoice_2.5kHz_EM_EM_0dB <- data_nochoice_2.5kHz_EM_EM %>% filter(TMR == "0dB")
      # Difference between signal alone and EM Control
      wilcox.test(data_nochoice_2.5kHz_EM_SA_0dB$Threshold, data_nochoice_2.5kHz_EM_EM_0dB$Threshold, paired = TRUE) 
      wilcox.test(data_nochoice_2.5kHz_EM_SA_0dB$Latency, data_nochoice_2.5kHz_EM_EM_0dB$Latency, paired = TRUE) 
      
      # Difference between IM Test and EM Control
      wilcox.test(data_nochoice_2.5kHz_IM_IM_0dB$Threshold, data_nochoice_2.5kHz_EM_EM_0dB$Threshold, alternative = "two.sided")
      wilcox.test(data_nochoice_2.5kHz_IM_IM_0dB$Latency, data_nochoice_2.5kHz_EM_EM_0dB$Latency, alternative = "two.sided")