---
title: "Phyllostomus discolor"
output: word_document

---

## Session Info:
```{r info}
# R version 3.6.2 (2019-12-12)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 7 x64 (build 7601) Service Pack 1
# 
# Matrix products: default
# 
# Random number generation:
#  RNG:     Mersenne-Twister 
#  Normal:  Inversion 
#  Sample:  Rounding 
#  
# locale:
# [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
# [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
# [5] LC_TIME=English_United States.1252    
# 
# attached base packages:
# [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#  [1] stringr_1.4.0   plyr_1.8.5      officer_0.3.6   flextable_0.5.8 magrittr_1.5   
#  [6] gridExtra_2.3   cowplot_1.0.0   dplyr_0.8.4     psyphy_0.1-9    ggstance_0.3.3 
# [11] ggplot2_3.2.1  
# 
# loaded via a namespace (and not attached):
#  [1] rstan_2.19.3       tidyselect_1.0.0   xfun_0.12          purrr_0.3.3       
#  [5] vctrs_0.2.3        colorspace_1.4-1   htmltools_0.4.0    stats4_3.6.2      
#  [9] loo_2.2.0          yaml_2.2.1         base64enc_0.1-3    utf8_1.1.4        
# [13] rlang_0.4.4        pkgbuild_1.0.6     pillar_1.4.3       glue_1.3.1        
# [17] withr_2.1.2        gdtools_0.2.1      uuid_0.1-2         matrixStats_0.55.0
# [21] lifecycle_0.1.0    munsell_0.5.0      gtable_0.3.0       zip_2.0.4         
# [25] evaluate_0.14      labeling_0.3       inline_0.3.15      knitr_1.28        
# [29] callr_3.4.2        ps_1.3.2           parallel_3.6.2     fansi_0.4.1       
# [33] Rcpp_1.0.3         scales_1.1.0       StanHeaders_2.19.0 farver_2.0.3      
# [37] systemfonts_0.1.1  packrat_0.5.0      digest_0.6.25      stringi_1.4.6     
# [41] processx_3.4.2     cli_2.0.1          tools_3.6.2        lazyeval_0.2.2    
# [45] tibble_2.1.3       crayon_1.3.4       pkgconfig_2.0.3    rsconnect_0.8.16  
# [49] xml2_1.2.2         data.table_1.12.8  prettyunits_1.1.1  assertthat_0.2.1  
# [53] rmarkdown_2.1      rstudioapi_0.11    R6_2.4.1           compiler_3.6.2  
```

```{r functions, include=FALSE, cache=F}
library(ggplot2)
library(ggstance)
library(psyphy)
library(dplyr)
library(cowplot)
library(gridExtra)
library(grid)
library(magrittr)
library(flextable)
# library(worded)
library(officer)
ref <- function(useName) {
require(stringr)
if(!exists(".refctr")) .refctr <- c(`_` = 0)
if(any(names(.refctr)==useName)) return(.refctr[useName])
type=str_split(useName,":")[[1]][1]
nObj <- sum(str_detect(names(.refctr),type))
useNum <- nObj + 1
newrefctr <- c(.refctr, useNum)
names(newrefctr)[length(.refctr) + 1] <- useName
assign(".refctr", newrefctr, envir=.rmdenvir)
return(useNum)
}
.rmdenvir = environment()
.refctr <- c(`_` = 0)
tab <- function(useNameT) {
require(stringr)
if(!exists(".refctrT")) .refctrT <- c(`_` = 0)
if(any(names(.refctrT)==useNameT)) return(.refctrT[useNameT])
type=str_split(useNameT,":")[[1]][1]
nObjT <- sum(str_detect(names(.refctrT),type))
useNumT <- nObjT + 1
newrefctrT <- c(.refctrT, useNumT)
names(newrefctrT)[length(.refctrT) + 1] <- useNameT
assign(".refctrT", newrefctrT, envir=.rmdenvir)
return(useNumT)
}
.rmdenvir = environment()
.refctr <- c(`_` = 0) 
```
```{r setup, include=FALSE, cache=F}
knitr::opts_chunk$set(fig.height=10,fig.width=14,echo = FALSE,warning=F,message=F,cache=T)
options(digits=3)

DE<-read.csv("Pdisc_Choice_Test.csv")
DE$Noise<-as.factor(DE$Noise)
DE$Noise<-relevel(DE$Noise,"Silent (Control)")
DE$MD<-DE$MD*2 # change from modulation depth to peak-to-peak measurement

DE$batID<-DE$ID
DE$batID<-as.character(DE$batID)
DE$batID[DE$batID=="Cartman"]<-"Bat A"
DE$batID[DE$batID=="Chiarah"]<-"Bat B"
DE$batID[DE$batID=="Anshin"]<-"Bat C"
DE$batID[DE$batID=="Mucki"]<-"Bat D"

```

```{r CorrModels,cache=T}

Cor1<-glm(Correct ~ Noise * MD +scale(Start),
               data=DE[DE$ID=="Cartman",], family=binomial(link="logit"))

Cor2<-glm(Correct ~ Noise * MD +scale(Start),
               data=DE[DE$ID=="Chiarah",], family=binomial(link="logit"))

Cor3<-glm(Correct ~ Noise * MD +scale(Start),
               data=DE[DE$ID=="Anshin",], family=binomial(link="logit"))

Cor4<-glm(Correct ~ Noise * MD +scale(Start),
               data=DE[DE$ID=="Mucki",], family=binomial(link="logit"))

CorrectTable<-data.frame("Bat"="A","Variable"=rownames(summary(Cor1)$coefficient),round(summary(Cor1)$coefficient,3))
letter<-c("A","B","C","D")
for(i in 2:4){
addTable<-data.frame("Bat"=letter[i],"Variable"=rownames(summary(get(paste0("Cor",i)))$coefficient),round(summary(get(paste0("Cor",i)))$coefficient,3))
CorrectTable<-rbind(CorrectTable,addTable)
}
rownames(CorrectTable)<-c(1:length(CorrectTable$Bat))
names(CorrectTable)[4:6]<-c("SE","Z value","p value")

CorrectTable$"p value"<-round(1 - (1-CorrectTable$"p value")^4,3) # bonferroni correction
CorrectTable$"p value"<-as.character(CorrectTable$"p value")
CorrectTable$"p value"<-ifelse(CorrectTable$"p value"=="0","<0.001",CorrectTable$"p value")
CorrectTable$Variable<-as.character(CorrectTable$Variable)
CorrectTable$Variable<-gsub(pattern = "Noise",replacement="",CorrectTable$Variable)
```
```{r CorrectScript, echo=F, cache=T}
Prop<- DE %>% dplyr::group_by(ID, Noise, MD) %>% 
  dplyr::summarize(
    N=n(),
    prop=sum(Correct==1)/dplyr::n(),
    sum=sum(Correct==1),
    sumInc=sum(Correct==0)
  )

bat<-as.character(unique(Prop$ID))[c(3,2,1,4)]
ID<-c("Bat A","Bat B","Bat C","Bat D")

i=1
for(i in 1:4){
set.seed(502)
  dat<-Prop[Prop$ID==bat[i],]

lMD<-unique(log(dat$MD))
dat$lMD<-log(dat$MD)


fitting <- function(df) {
  model <- glm((df$sum/df$N) ~ lMD, family = binomial(mafc.probit(2)), weights=df$N)
  xseq <- seq(min(lMD), max(lMD), len = 1000)  #I used, for example, a 1000 points
  yseq <- predict(model, data.frame(lMD=xseq), type="response")
  data.frame(xseq, yseq)
}
library(plyr)
curves <-ddply(dat, .(Noise), fitting)

One<-ifelse(i==2, T, F)
library(ggplot2)
p1 <- ggplot() + 
  geom_point(data = dat, aes(x = exp(lMD), y = prop, color=Noise), show.legend = One, size=4)+ 
  coord_cartesian(ylim=c(0.4,1), xlim=c(2,32))+
  theme_classic()+
  geom_hline(yintercept=0.5,linetype="dashed")+
  geom_hline(yintercept=1,linetype="dashed")+
  scale_x_log10(breaks=c(sort(unique(DE$MD))))+
  scale_y_continuous(breaks=seq(.4,1,by=.10),
                       labels=seq(40,100,by=10))+
  theme(axis.text=element_text(size=22),
        legend.text = element_text(size=20),
        plot.title = element_text(size=26))+
  labs(title=ID[i])


## Legends ####
if(i==1){
          p1<- p1+theme(axis.title=element_blank(),
                    axis.text.x=element_blank(),
                         axis.ticks.length.y = unit(.25,"cm"),
                          axis.ticks.y = element_line(size=2),
                         panel.grid.major.y = element_line(colour = "grey", size = .5)) 
} else if(i==2) {
          p1<- p1+theme(axis.title=element_blank(),
                legend.position = c(.78,.3),
                legend.title = element_blank(),
                axis.text=element_blank(),
                panel.grid.major.y = element_line(colour = "grey", size = .5)) 
} else if(i==4) {
          p1<- p1+theme(axis.title=element_blank(),
                axis.text.y=element_blank(),
                 axis.text.x = element_text(angle=45,vjust=0.6),
                         axis.ticks.length.x = unit(.25,"cm"),
                          axis.ticks.x = element_line(size=2),
                         panel.grid.major.y = element_line(colour = "grey", size = .5)) 
} else {  
           p1<- p1+theme(axis.title=element_blank(),
                          axis.text.x = element_text(angle=45,vjust=0.6),
                         axis.ticks.length = unit(.25,"cm"),
                          axis.ticks = element_line(size=2),
                         panel.grid.major.y = element_line(colour = "grey", size = .5)) 
  }

## Add curves ####
p1<-p1 + geom_line(data = curves, aes(x = exp(xseq), y = yseq, color = Noise), size=1.2, alpha=0.5, show.legend = One)

threshold<-function(df){
model <- glm((df$sum/df$N) ~ lMD, family = binomial(mafc.probit(2)), weights=df$N)
(m <- -coef(model)[[1]]/coef(model)[[2]])
(std <- 1/coef(model)[[2]])
p <- 0.7
(thre <- qnorm((p - 0.5)/0.5, m, std))
threshold <- data.frame(p, thre)
}
thresholds <- ddply(dat, .(Noise), threshold)

if(i==1){
  Threshold_Table<-thresholds
} else {
  Threshold_Table<-rbind(Threshold_Table,thresholds)
}

p1<-p1 + geom_linerange(data = thresholds,size=.8,show.legend = F,
                        aes(x = exp(thre), ymin = 0, ymax = p,  color = Noise),alpha=0.5)


## error bar
library("plyr")

sampling <- function(x) {
  funct <- function(df) {
    n <- length(df$lMD)
    size <- df$sum + df$sumInc
    lMD <- df$lMD
    nCorrect <- rbinom(n, size, prob = df$prop)
    nIncorrect <- size - nCorrect
    data.frame(lMD, nCorrect, nIncorrect)
  }
  ddply(dat, .(Noise), funct)
}
fake.dat <- ddply(data.frame(sample = 1:1000), .(sample), sampling)


calculate.thre <- function(df) {
  tryCatch({
    model <- glm(cbind(df$nCorrect,df$nIncorrect) ~ lMD, family = binomial(mafc.probit(2)))
    m <- -coef(model)[[1]]/(coef(model)[[2]])
    std <- 1/(coef(model)[[2]])
    p <- 0.7
    thre <- qnorm((p - 0.5)/0.5, m, std)
    data.frame(thre)
  }, error = function(e) message(paste("Fitting error in sample: ", unique(df$sample), 
                                       "\n")), warning = function(w) message(paste("Fitting warning in sample: ", 
                                                                                   unique(df$sample), "\n")))
}
fake.thresholds <- ddply(fake.dat, .(sample, Noise), calculate.thre)
p <- 0.7
calculate.ci <- function(df) { # function to grab the 95% of the middle data
  conf.int <- quantile(df$thre, c(0.025, 0.975), na.rm = T) # grab data that are at 0.025 and 0.975 percentiles
  inf <- conf.int[[1]] # assign lower bound of 95% CI to object 
  sup <- conf.int[[2]] # assign upper bound of 95% CI to object 
  data.frame(inf, sup)
}
cis <- ddply(fake.thresholds, .(Noise), calculate.ci)

if(i==1){
  CI_thresh<-cis
} else {
  CI_thresh<-rbind(CI_thresh,cis)
}


p2<-p1 +  geom_errorbarh(data = cis, aes(y = p+c(0.0075,-0.00125,-0.0075,0.00125), xmax = exp(sup), xmin = exp(inf), color=Noise),alpha=0.7,height = 0.03, size=1,show.legend = F)

assign(noquote(bat[i]),p2)
}

Threshold_Table$ID<-c(rep("A",4),rep("B",4),rep("C",4),rep("D",4))
Sigtable<-cbind(Threshold_Table,CI_thresh)
Sigtable<-Sigtable[,-5]
Sigtable$thre<-round(exp(Sigtable$thre),2)
Sigtable$inf<-round(exp(Sigtable$inf),2)
Sigtable$sup<-round(exp(Sigtable$sup),2)

Sigtable<-Sigtable[order(Sigtable[,1]),]

names(Sigtable)<-c("Noise","p","Threshold","ID","Lower 95% CI","Upper 95% CI")

  Performance<-plot_grid(Chiarah,Cartman,Anshin,Mucki, align="vh")

y.grob <- textGrob("Correct (%)\n", 
                   gp=gpar(fontface="bold", fontsize=26), rot=90)

x.grob <- textGrob("\nRipple height (mm)", 
                   gp=gpar(fontface="bold", fontsize=26))

```
```{r TimeMeans,cache=T}
MeanTime<-DE %>% group_by(batID,Noise) %>%  dplyr::summarise(
  Mean=round(mean(Time,na.rm=T),2),
  SD=sd(Time,na.rm=T)
)
MeanTime<-as.data.frame(MeanTime)
```
```{r ModelTime,cache=T}
name<-sort(as.character(unique(DE$batID)))
for(i in 1:4){
  glm<-glm(log(Time) ~ Noise * MD +scale(Start),
           data=DE[DE$batID==name[i],], 
           family=inverse.gaussian(link="identity"))
  assign(paste0("glm",i),glm)
}
TimeTable<-data.frame("Bat"="A","Variable"=rownames(summary(glm1)$coefficient),round(summary(glm1)$coefficient,3))
letter<-c("A","B","C","D")
for(i in 2:4){
addTable<-data.frame("Bat"=letter[i],"Variable"=rownames(summary(get(paste0("glm",i)))$coefficient),round(summary(get(paste0("glm",i)))$coefficient,3))
TimeTable<-rbind(TimeTable,addTable)
}

rownames(TimeTable)<-c(1:length(TimeTable$Bat))
names(TimeTable)[4:6]<-c("SE","Z value","p value")

TimeTable$"p value"<-round(1 - (1-TimeTable$"p value")^4,3) # bonferroni correction
TimeTable$"p value"<-as.character(TimeTable$"p value")
TimeTable$"p value"<-ifelse(TimeTable$"p value"=="0","<0.001",TimeTable$"p value")
TimeTable$Variable<-as.character(TimeTable$Variable)
TimeTable$Variable<-gsub(pattern = "Noise",replacement="",TimeTable$Variable)
```
```{r ModelAbort,cache=T}
ID<-c("Bat A","Bat B","Bat C","Bat D")
for(i in 1:4){
  glm<-glm(abort ~ Noise+MD +scale(Start),
           data=DE[DE$batID==ID[i],], 
           family=binomial(link = "logit"))
  assign(paste0("glm",i),glm)
}

AbortTable<-data.frame("Bat"="A","Variable"=rownames(summary(glm1)$coefficient),round(summary(glm1)$coefficient,3))
letter<-c("A","B","C","D")
for(i in 2:4){
  addTable<-data.frame("Bat"=letter[i],"Variable"=rownames(summary(get(paste0("glm",i)))$coefficient),round(summary(get(paste0("glm",i)))$coefficient,3))
  AbortTable<-rbind(AbortTable,addTable)
}

rownames(AbortTable)<-c(1:length(AbortTable$Bat))
names(AbortTable)[4:6]<-c("SE","Z value","p value")

AbortTable$"p value"<-round(1 - (1-AbortTable$"p value")^4,3) # bonferroni correction
AbortTable$"p value"<-as.character(AbortTable$"p value")
AbortTable$"p value"<-ifelse(AbortTable$"p value"=="0","<0.001",AbortTable$"p value")
AbortTable$Variable<-as.character(AbortTable$Variable)
AbortTable$Variable<-gsub(pattern = "Noise",replacement="",AbortTable$Variable)
```
```{r AbortProp,cache=T}
Prop<- 
DE %>% dplyr::group_by(batID, Noise) %>% 
  dplyr::summarize(
    N=n(),
    AbortProp=sum(abort==1)/(N+sum(abort==1)),
    Aborts=sum(abort==1),
    sd=sd(abort),
    se=sd/sqrt(N),
    CI=1.96*se
  )

```
```{r CAsetup,cache=T}
CA<-read.csv("Pdisc_CA.csv")

CA$mod_depth<-2*CA$mod_depth # change to peak-to-peak measurement
CAc<-CA %>% group_by(tot_trial_c,mod_depth,bat_num,
                     Treat, start_t,corr,latency) %>% 
  dplyr::summarise(
    count=n(),
    dur=mean(dur..ms.),
    RMS=mean(RMS),
    Peak=mean(peakF_spg..kHz.),
    Centroid=mean(Fcentroid..kHz.)
  )

CAc$Treat<-as.factor(CAc$Treat)
CAc$Treat<-relevel(CAc$Treat,"Silent (Control)")
levels(CAc$Treat)<-c("Silence","Non-overlap","Smooth-overlap","Sparse-overlap")
```
```{r durMod,cache=T}
bat<-sort(unique(CAc$bat_num))
letter<-c("A","B","C","D")
for(i in 1:4){
assign(paste0("dur",letter[i]), 
       lm(log(dur) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],])
)
assign(paste0("DurCoef",letter[i]), 
       summary(lm(log(dur) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],]))$coefficients
       )
}
  
DurTable<-data.frame("Bat"="A","Variable"=rownames(DurCoefA),round(DurCoefA,3))
letter<-c("A","B","C","D")
for(i in 2:4){
addTable<-data.frame("Bat"=letter[i],
                  "Variable"=rownames(get(paste0("DurCoef",letter[i]))),
                    round(get(paste0("DurCoef",letter[i])),3))
DurTable<-rbind(DurTable,addTable)
}

rownames(DurTable)<-c(1:length(DurTable$Bat))
names(DurTable)[4:6]<-c("SE","Z value","p value")

DurTable$"p value"<-round(1 - (1-DurTable$"p value")^4,3) # bonferroni correction
DurTable$"p value"<-as.character(DurTable$"p value")
DurTable$"p value"<-ifelse(DurTable$"p value"=="0","<0.001",DurTable$"p value")
DurTable$Variable<-as.character(DurTable$Variable)
DurTable$Variable<-gsub(pattern ="Noise",replacement="",DurTable$Variable)
DurTable$Variable<-gsub(pattern ="Treat",replacement="",DurTable$Variable)

```
```{r RMSMod,cache=T}
bat<-sort(unique(CAc$bat_num))
letter<-c("A","B","C","D")
for(i in 1:4){
  assign(paste0("RMS",letter[i]), 
         lm(log(abs(RMS)) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],])
  )
  assign(paste0("RMSCoef",letter[i]), 
         summary(lm(log(abs(RMS)) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],]))$coefficients
  )
}

RMSTable<-data.frame("Bat"="A","Variable"=rownames(RMSCoefA),round(RMSCoefA,3))
letter<-c("A","B","C","D")
for(i in 2:4){
  addTable<-data.frame("Bat"=letter[i],
                       "Variable"=rownames(get(paste0("RMSCoef",letter[i]))),
                       round(get(paste0("RMSCoef",letter[i])),3))
  RMSTable<-rbind(RMSTable,addTable)
}

rownames(RMSTable)<-c(1:length(RMSTable$Bat))
names(RMSTable)[4:6]<-c("SE","Z value","p value")

RMSTable$"p value"<-round(1 - (1-RMSTable$"p value")^4,3) # bonferroni correction
RMSTable$"p value"<-as.character(RMSTable$"p value")
RMSTable$"p value"<-ifelse(RMSTable$"p value"=="0","<0.001",RMSTable$"p value")
RMSTable$Variable<-as.character(RMSTable$Variable)
RMSTable$Variable<-gsub(pattern ="Noise",replacement="",RMSTable$Variable)
RMSTable$Variable<-gsub(pattern ="Treat",replacement="",RMSTable$Variable)
```
```{r ModPeak,cache=T}
bat<-sort(unique(CAc$bat_num))
letter<-c("A","B","C","D")
for(i in 1:4){
  assign(paste0("Peak",letter[i]), 
         lm(log(Peak) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],])
  )
  assign(paste0("PeakCoef",letter[i]), 
         summary(lm(log(Peak) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],]))$coefficients
  )
}


PeakTable<-data.frame("Bat"="A","Variable"=rownames(PeakCoefA),round(PeakCoefA,3))
letter<-c("A","B","C","D")
for(i in 2:4){
  addTable<-data.frame("Bat"=letter[i],
                       "Variable"=rownames(get(paste0("PeakCoef",letter[i]))),
                       round(get(paste0("PeakCoef",letter[i])),3))
  PeakTable<-rbind(PeakTable,addTable)
}

rownames(PeakTable)<-c(1:length(PeakTable$Bat))
names(PeakTable)[4:6]<-c("SE","Z value","p value")

PeakTable$"p value"<-round(1 - (1-PeakTable$"p value")^4,3) # bonferroni correction
PeakTable$"p value"<-as.character(PeakTable$"p value")
PeakTable$"p value"<-ifelse(PeakTable$"p value"=="0","<0.001",PeakTable$"p value")
PeakTable$Variable<-as.character(PeakTable$Variable)
PeakTable$Variable<-gsub(pattern ="Noise",replacement="",PeakTable$Variable)
PeakTable$Variable<-gsub(pattern ="Treat",replacement="",PeakTable$Variable)
```
```{r ModCent,cache=T}
bat<-sort(unique(CAc$bat_num))
letter<-c("A","B","C","D")
for(i in 1:4){
  assign(paste0("Cent",letter[i]), 
         lm(log(Centroid) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],])
  )
  assign(paste0("CentCoef",letter[i]), 
         summary(lm(log(Centroid) ~ Treat*mod_depth+scale(start_t),data=CAc[CAc$bat_num==bat[i],]))$coefficients
  )
}

CentTable<-data.frame("Bat"="A","Variable"=rownames(CentCoefA),round(CentCoefA,3))
letter<-c("A","B","C","D")
for(i in 2:4){
  addTable<-data.frame("Bat"=letter[i],
                       "Variable"=rownames(get(paste0("CentCoef",letter[i]))),
                       round(get(paste0("CentCoef",letter[i])),3))
  CentTable<-rbind(CentTable,addTable)
}

rownames(CentTable)<-c(1:length(CentTable$Bat))
names(CentTable)[4:6]<-c("SE","Z value","p value")

CentTable$"p value"<-round(1 - (1-CentTable$"p value")^4,3) # bonferroni correction
CentTable$"p value"<-as.character(CentTable$"p value")
CentTable$"p value"<-ifelse(CentTable$"p value"=="0","<0.001",CentTable$"p value")
CentTable$Variable<-as.character(CentTable$Variable)
CentTable$Variable<-gsub(pattern ="Noise",replacement="",CentTable$Variable)
CentTable$Variable<-gsub(pattern ="Treat",replacement="",CentTable$Variable)
```

## Results

All four bats learned to discriminate the smooth disc from the rippled disc with the highest ripples in silence (88-100% correct), and showed reduced discrimination performance with decreasing ripple height (**`r paste("Figure", ref("fig: FigCorr"))`**, orange line; logistic model p `r CorrectTable[CorrectTable$Bat=="A",][5,6]`; **`r paste("Table", tab("fig: TableCorr"))`**). In silence, performance dropped below our 70% threshold level for ripple heights between `r round(min(Sigtable$Threshold[1:4]),1)` mm - `r round(max(Sigtable$Threshold[1:4]),1)` mm (**`r paste("Figure", ref("fig: FigCorr"))`**, **`r paste("Table", tab("fig: TableSigmoid"))`**, see thresholds). Under noise treatments, the discrimination performance of bats A and B (hereafter 'coping' bats) did not change (**`r paste("Table", tab("fig: TableCorr"))`**), and the 95% confidence intervals for the thresholds in noise overlapped with those in silence. In contrast, discrimination performance decreased for bats C and D (hereafer 'non-coping' bats) in both smooth-overlapping (z = `r round(CorrectTable[CorrectTable$Bat=="C",][8,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="C",][8,6]`; z = `r round(CorrectTable[CorrectTable$Bat=="D",][8,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="D",][8,6]`
) and sparse-overlapping (z = `r round(CorrectTable[CorrectTable$Bat=="C",][9,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="C",][9,6]`; z = `r round(CorrectTable[CorrectTable$Bat=="D",][9,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="D",][9,6]`
) noise treatments (blue and purple lines in **`r paste("Figure", ref("fig: FigCorr"))`** respectively). The same is true for the non-overlapping noise treatment for bat C (z = `r round(CorrectTable[CorrectTable$Bat=="C",][7,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="C",][7,6]`), yet not for bat D (z = `r round(CorrectTable[CorrectTable$Bat=="D",][7,5],1)`, p = `r  CorrectTable[CorrectTable$Bat=="D",][7,6]`).

Time to complete trials differed between noise treatments (**`r paste("Figure", ref("fig: FigureTime"))`**). However, responses differed among bats. Both bats A and D were quicker to make a decision during smooth-overlapping noise compared to silence (mean times of bat A and D under noise and silence, respectively: `r round(MeanTime[MeanTime$batID=="Bat A",][3,3],1)` s vs. `r round(MeanTime[MeanTime$batID=="Bat A",][1,3],1)` s (A) and `r round(MeanTime[MeanTime$batID=="Bat D",][3,3],1)` s vs. `r round(MeanTime[MeanTime$batID=="Bat D",][1,3],1)` s (D); z = `r round(TimeTable[TimeTable$Bat=="A",][3,5],1)`, p = `r TimeTable[TimeTable$Bat=="A",][3,6]`; z = `r round(TimeTable[TimeTable$Bat=="D",][3,5],1)`, p = `r TimeTable[TimeTable$Bat=="D",][3,6]`, bats A and D respectively; **`r paste("Table", tab("fig: TableTime"))`**). However, bat C took longer to complete trials during sparse-overlapping noise (`r round(MeanTime[MeanTime$batID=="Bat C",][4,3],1)` s vs. `r round(MeanTime[MeanTime$batID=="Bat C",][1,3],1)` s; z = `r round(TimeTable[TimeTable$Bat=="C",][4,5],1)`, p = `r TimeTable[TimeTable$Bat=="C",][4,6]`), while noise treatments did not affect the trial time of bat B (**`r paste("Table", tab("fig: TableTime"))`**).


There were a total of `r sum(Prop$Aborts)` aborted trials (bats A: `r sum(Prop$Aborts[Prop$batID=="Bat A"])`; B: `r sum(Prop$Aborts[Prop$batID=="Bat B"])`; C: `r sum(Prop$Aborts[Prop$batID=="Bat C"])`; D: `r sum(Prop$Aborts[Prop$batID=="Bat D"])`) of `r sum(Prop$Aborts)+length(DE$Correct)` total trials (`r round(sum(Prop$Aborts)/(sum(Prop$Aborts)+length(DE$Correct))*100,1)` %). Compared to silence, both bats B and C significantly aborted more trials under both non-overlapping (z = `r format(round(AbortTable[AbortTable$Bat=="B",][2,5],1), nsmall = 1)`, p = `r AbortTable[AbortTable$Bat=="B",][2,6]`; z = `r format(round(AbortTable[AbortTable$Bat=="C",][2,5],1), nsmall = 1)`, p = `r AbortTable[AbortTable$Bat=="C",][2,6]`) and smooth-overlapping noise (z = `r round(AbortTable[AbortTable$Bat=="B",][3,5],1)`, p = `r AbortTable[AbortTable$Bat=="B",][3,6]`; z = `r round(AbortTable[AbortTable$Bat=="C",][3,5],1)`, p = `r AbortTable[AbortTable$Bat=="C",][3,6]`). Bats B and D aborted more trials under sparse-overlapping noise compared to silence (z = `r round(AbortTable[AbortTable$Bat=="B",][4,5],1)`, p = `r AbortTable[AbortTable$Bat=="B",][4,6]`; z = `r round(AbortTable[AbortTable$Bat=="D",][4,5],1)`, p = `r AbortTable[AbortTable$Bat=="D",][4,6]`; **`r paste("Figure", ref("fig: FigAbort"))`**; **`r paste("Table", tab("fig: TableAbort"))`**). 


In order to understand how bats may be altering echolocation calls to cope with noise, we recorded `r length(CA$NoiseRMS)` individual echolocation calls across the `r length(DE$Correct)` trials. Call characteristics were measured and averaged by trial for analysis. 

The mean call duration was `r round(mean(CA$dur),2)` ms (bats A: `r round(mean(CA$dur[CA$bat_num=="bat A"]),2)` ms; B: `r round(mean(CA$dur[CA$bat_num=="bat B"]),2)` ms; C: `r round(mean(CA$dur[CA$bat_num=="bat C"]),2)` ms; D: `r round(mean(CA$dur[CA$bat_num=="bat D"]),2)` ms). All bats increased call duration under smooth-overlapping noise, however effect sizes varied. Coping bats increased call duration by `r exp(DurCoefA[1,1]+DurCoefA[3,1])-exp(DurCoefA[1,1])` ms and `r exp(DurCoefB[1,1]+DurCoefB[3,1])-exp(DurCoefB[1,1])` ms, while non-coping bats only increased call duration by `r exp(DurCoefC[1,1]+DurCoefC[3,1])-exp(DurCoefC[1,1])` ms and `r exp(DurCoefD[1,1]+DurCoefD[3,1])-exp(DurCoefD[1,1])` ms (t = `r format(round(DurCoefA[3,3],1),nsmall=1)`, p = `r DurTable[DurTable$Bat=="A",][3,6]`; t = `r round(DurCoefB[3,3],1)`, p = `r DurTable[DurTable$Bat=="B",][3,6]`; t = `r round(DurCoefC[3,3],1)`, p = `r DurTable[DurTable$Bat=="C",][3,6]`; t = `r round(DurCoefD[3,3],1)`, p = `r DurTable[DurTable$Bat=="D",][3,6]`; bats A, B, C, D respectively; **`r paste("Table", tab("fig: TableDur"))`**). Similarly, coping bats increased call duration in sparse-overlapping noise (increase of `r exp(DurCoefA[1,1]+DurCoefA[4,1])-exp(DurCoefA[1,1])` ms and `r exp(DurCoefB[1,1]+DurCoefB[4,1])-exp(DurCoefB[1,1])` ms, A and B respectively), while non-coping bats did not (t = `r round(DurCoefA[4,3],1)`, p = `r DurTable[DurTable$Bat=="A",][4,6]`; t = `r round(DurCoefB[4,3],1)`, p = `r DurTable[DurTable$Bat=="B",][4,6]`; t = `r round(DurCoefC[4,3],1)`, p = `r DurCoefC[4,4]`; t = `r round(DurCoefD[4,3],1)`, p = `r DurCoefD[4,4]`; bats A, B, C, D respectively; **`r paste("Figure", ref("fig: FigDur"))`**). Oddly, bats B and C decreased call duration by `r exp(DurCoefB[1,1]+DurCoefB[2,1])-exp(DurCoefB[1,1])` ms and `r exp(DurCoefC[1,1]+DurCoefC[2,1])-exp(DurCoefC[1,1])` ms in non-overlapping noise relative to silence (t = `r round(DurCoefB[2,3],1)`, p = `r DurTable[DurTable$Bat=="B",][2,6]`; t = `r round(DurCoefC[2,3],1)`, p = `r DurTable[DurTable$Bat=="C",][2,6]`; bats B and C respectively).

Relative to silence, all bats increased sonar amplitude (RMS) in smooth-overlapping noise (bats A: `r abs(round(exp(RMSCoefA[1,1]+RMSCoefA[3,1])-exp(RMSCoefA[1,1]),1))` dB; B:  `r abs(round(exp(RMSCoefB[1,1]+RMSCoefB[3,1])-exp(RMSCoefB[1,1]),1))` dB; C: `r abs(round(exp(RMSCoefC[1,1]+RMSCoefC[3,1])-exp(RMSCoefC[1,1]),1))` dB; D: `r abs(round(exp(RMSCoefD[1,1]+RMSCoefD[3,1])-exp(RMSCoefD[1,1]),1))` dB) and in sparse-overlapping noise (bats A: `r abs(round(exp(RMSCoefA[1,1]+RMSCoefA[4,1])-exp(RMSCoefA[1,1]),1))` dB; B: `r abs(round(exp(RMSCoefB[1,1]+RMSCoefB[4,1])-exp(RMSCoefB[1,1]),1))` dB; C: `r abs(round(exp(RMSCoefC[1,1]+RMSCoefC[4,1])-exp(RMSCoefC[1,1]),1))` dB; D: `r abs(round(exp(RMSCoefD[1,1]+RMSCoefD[4,1])-exp(RMSCoefD[1,1]),1))` dB; **`r paste("Figure", ref("fig: FigRMS"))`**; **`r paste("Table", tab("fig: TableRMS"))`**; results for bats A, B, C, D respectively in smooth-overlapping: t = `r round(RMSCoefA[3,3],1)`; t = `r round(RMSCoefB[3,3],1)`; t = `r round(RMSCoefC[3,3],1)`; t = `r round(RMSCoefD[3,3],1)`; and in sparse-overlapping: t = `r round(RMSCoefA[4,3],1)`; t = `r round(RMSCoefB[4,3],1)`; t = `r round(RMSCoefC[4,3],1)`; t = `r round(RMSCoefD[4,3],1)`; p for all tests < 0.001).  Bat A also had a significant, but much smaller, increase in relative amplitude during the non-overlapping noise treatment (increase of `r abs(round(exp(RMSCoefA[1,1]+RMSCoefA[2,1])-exp(RMSCoefA[1,1]),1))` dB, t = `r round(RMSCoefA[2,3],1)`, p = `r RMSTable[RMSTable$Bat=="A",][2,6]`).

The mean peak frequency was `r round(mean(CAc$Peak),1)` kHz (bats A: `r round(mean(CAc$Peak[CAc$bat_num=="bat A"]),1)` kHz; B: `r round(mean(CAc$Peak[CAc$bat_num=="bat B"]),1)` kHz; C: `r round(mean(CAc$Peak[CAc$bat_num=="bat C"]),1)` kHz; D: `r round(mean(CAc$Peak[CAc$bat_num=="bat D"]),1)` kHz). Bat A increased peak frequency in non-overlapping noise by `r round(exp(PeakCoefA[1,1]+PeakCoefA[2,1])-exp(PeakCoefA[1,1]),1)` kHz, and decreased peak frequency in overlapping noise by `r round(exp(PeakCoefA[1,1]+PeakCoefA[3,1])-exp(PeakCoefA[1,1]),1)` kHz (t = `r round(PeakCoefA[2,3],1)`, p = `r PeakTable[PeakTable$Bat=="A",][2,6]`; t = `r round(PeakCoefA[3,3],1)`, p = `r PeakTable[PeakTable$Bat=="A",][3,6]`). Bat C increased peak frequency by `r round(exp(PeakCoefC[1,1]+PeakCoefC[4,1])-exp(PeakCoefC[1,1]),1)` kHz only in sparse-overlapping noise (t = `r round(PeakCoefC[4,3],1)`, p = `r PeakTable[PeakTable$Bat=="C",][4,6]`; **`r paste("Figure", ref("fig: FigPeak"))`**). Bats B and D did not change the peak frequency in any of the noise treatments (**`r paste("Table", tab("fig: TablePeak"))`**).

The mean frequency centroid was `r round(mean(CAc$Centroid),1)` kHz (bats A: `r round(mean(CAc$Centroid[CAc$bat_num=="bat A"]),1)` kHz; B: `r round(mean(CAc$Centroid[CAc$bat_num=="bat B"]),1)` kHz; C: `r round(mean(CAc$Centroid[CAc$bat_num=="bat C"]),1)` kHz; D: `r round(mean(CAc$Centroid[CAc$bat_num=="bat D"]),1)` kHz). Bats A, B, and C all decreased the frequency centroid in overlapping noise by `r format(round(exp(CentCoefA[1,1]+CentCoefA[3,1])-exp(CentCoefA[1,1]),1),nsmall=1)` kHz, `r round(exp(CentCoefB[1,1]+CentCoefB[3,1])-exp(CentCoefB[1,1]),1)` kHz, and `r round(exp(CentCoefC[1,1]+CentCoefC[3,1])-exp(CentCoefC[1,1]),1)` kHz respectively (A: t = `r round(CentCoefA[3,3],1)`, p = `r CentTable[CentTable$Bat=="A",][3,6]`; B:  t = `r round(CentCoefB[3,3],1)`, p = `r CentTable[CentTable$Bat=="B",][3,6]`; C: , t = `r round(CentCoefC[3,3],1)`, p = `r CentTable[CentTable$Bat=="C",][3,6]`), while only bat A decreased the centroid in sparse noise (`r round(exp(CentCoefA[1,1]+CentCoefA[4,1])-exp(CentCoefA[1,1]),1)` kHz, t = `r round(CentCoefA[4,3],1)`, p = `r CentTable[CentTable$Bat=="A",][4,6]`; **`r paste("Figure", ref("fig: FigCent"))`**; **`r paste("Table", tab("fig: TableCent"))`**). Bat A also increased the frequency centroid non-overlapping noise by `r round(exp(CentCoefA[1,1]+CentCoefA[2,1])-exp(CentCoefA[1,1]),1)` kHz, while non-coping bats decreased the frequency centroid in this treatment by `r round(exp(CentCoefC[1,1]+CentCoefC[2,1])-exp(CentCoefC[1,1]),1)` kHz and `r round(exp(CentCoefD[1,1]+CentCoefD[2,1])-exp(CentCoefD[1,1]),1)` kHz respectively (A: t = `r format(round(CentCoefA[2,3],1),nsmall=1)`, p = `r CentTable[CentTable$Bat=="A",][2,6]`; C: t = `r format(round(CentCoefC[2,3],1),nsmall=1)`, p = `r CentTable[CentTable$Bat=="C",][2,6]`; D: t = `r round(CentCoefD[2,3],1)`, p = `r CentTable[CentTable$Bat=="D",][2,6]`).


## Figures:

```{r FigCorr,fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigCorr"),"Discrimination performance of four bats as a function of disc ripple height during silent control and three noise treatments.** Asterisks denote that the interaction between noise treatment and modulation depth (rel. silent controls; orange lines) differs significantly in generalized linear models. Discrimination performance of 'coping' bats (A and B) did not differ between silence and noise treatments, whereas discrimination performance of 'non-coping' bats (C and D) was reduced under both overlapping noise treatments (blue and purple lines), while discrimination performance decreased in non-overlapping noise (green lines) for bat C only.")}
grid.arrange(arrangeGrob(Performance,
                           left = y.grob, bottom = x.grob)   )

```

```{r FigureTime,fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigureTime"),"Trial duration by noise treatment.** ...")}
ggplot(data=DE, aes(y=Time,x=Noise,fill=Noise) )+
         geom_boxplot(width=0.5)+
         theme_minimal()+
         facet_wrap(~batID)+
         theme(legend.position = "none",
               axis.line=element_line(size=1),
               axis.text = element_text(size=22),
               axis.text.x = element_text(angle = 45,vjust = .6),
               axis.title = element_text(size=26),
               strip.background =element_rect(color="white", fill="white"),
               strip.text = element_text(size=26, hjust = 0.025),
               panel.grid.major = element_line(color="grey80")
               )+
        scale_x_discrete(labels=c("Control","Smooth non-OL","Smooth OL","Sparse OL"))+
         stat_summary(fun.y=median, geom="point", shape=23, size=2)+
         scale_y_log10()+
  labs(y="Time (s)\n", x="Noise treatment")
```

```{r FigAbort,fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigAbort"),"The proportion of aborted trials plotted by noise treatment.** Both bats B and C had a significant increase in the proportion of aborted trials for both non-overlapping and overlapping noise. Yet, bats B and D had a significant increase in the proportion of aborted trials for sparse overlapping noise.")}

ggplot( aes(x=Noise,y=AbortProp,fill=Noise),data=Prop)+
  facet_wrap(~batID)+
  geom_col()+
  geom_errorbar(aes(ymin=AbortProp-CI,ymax=AbortProp+CI),width=0.25)+
  theme_classic()+
  theme(legend.position = "none",
    axis.text = element_text(size=22),
    axis.ticks.length = unit(.25,"cm"),
    axis.ticks = element_line(size=2),
    axis.title = element_text(size=26),
    strip.background =element_rect(color="white", fill="white"),
    strip.text = element_text(size=26, hjust = 0.025), 
    panel.grid.major.y = element_line(colour = "grey", size = .5),
    axis.text.x = element_text(angle = 45,vjust = .6)
    )+
  scale_y_continuous(breaks=c(0,.05,.10,.15,.20),labels=c(0,5,10,15,20))+
  coord_cartesian(ylim=c(0,0.23))+
  scale_x_discrete(labels=c("Control","Smooth non-OL","Smooth OL","Sparse OL"))+
  labs(y="Aborted trials (%)\n",x="\nNoise treatment")

```

```{r FigDur,fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigDur"),"Echolocation call duration.** ")}
CAc$bat_num<-as.character(CAc$bat_num)
CAc$bat_num<-gsub("bat","Bat",CAc$bat_num)
ggplot(aes(x=Treat,y=dur,fill=Treat),data=CAc)+
  geom_boxplot(width=0.5)+
  facet_wrap(~bat_num)+
  theme_minimal()+
  labs(x="Noise treatment",y="Call duration (ms)\n")+
  coord_cartesian(ylim=c(0.32,0.75))+
  theme(legend.position ="none",
        axis.line=element_line(size=1),
        panel.grid.major = element_line(color="grey80"),
        strip.background =element_rect(fill="white", color="white"),
        strip.text = element_text(size=26, hjust = 0.025),
        axis.text=element_text(size=22),
        axis.text.x=element_text(angle=45,vjust=0.55),
        axis.title=element_text(size=26)
  )
```

```{r FigRMS,fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigRMS"),"Echolocation call amplitude (RMS).** ")}
ggplot(aes(x=Treat,y=RMS,fill=Treat ),data=CAc)+
  geom_boxplot(width=.5)+ 
  facet_wrap(~bat_num)+
  theme_minimal()+
  labs(x="Noise treatment",y="Call RMS (dB FS)\n")+
    theme(legend.position = "none",
        axis.line=element_line(size=1),
        panel.grid.major = element_line(color="grey80"),
        strip.background = element_rect(fill="white", color="white"),
        strip.text = element_text(size=26, hjust = 0.025),
        axis.text=element_text(size=22),
        axis.text.x=element_text(angle=45,vjust=0.5),
        axis.title=element_text(size=26)
  )+
  coord_cartesian(ylim=c(-35,-5))

```

```{r FigPeak, fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigPeak"),"Echolocation signal peak frequency.** ")}
ggplot(aes(x=Treat,y=Peak,fill=Treat),data=CAc)+
  geom_boxplot(width=0.5)+
  facet_wrap(~bat_num)+ 
  theme_minimal()+
  labs(x="Noise treatment",y="Call peak frequency\n")+
  scale_y_log10()+
  coord_cartesian(ylim=c(55,90))+
  theme(legend.position = "none",
                axis.line=element_line(size=1),
        panel.grid.major = element_line(color="grey80"),
        strip.background =element_rect(fill="white",color="white"),
        strip.text = element_text(size=26, hjust = 0.025),
        axis.text=element_text(size=22),
        axis.text.x=element_text(angle=45,vjust=0.5),
        axis.title=element_text(size=26)
  )
```

```{r FigCent, fig.fullwidth=T,fig.cap=paste("**Figure",ref("fig: FigCent"),"Echolocation signal frequency centroid.** ")}
ggplot(aes(x=Treat,y=Centroid,fill=Treat), data=CAc)+
  geom_boxplot(width=0.5)+ 
  facet_wrap(~bat_num)+
  theme_minimal()+
  labs(x="Treatment",y="Call frequency Centroid\n")+
  scale_y_log10()+
  coord_cartesian(ylim=c(60,75))+
  guides(fill=guide_legend(ncol=2))+
  theme(legend.position = "none",
        axis.line=element_line(size=1),
        panel.grid.major = element_line(color="grey80"),
        strip.background =element_rect(fill="white",color="white"),
        strip.text = element_text(size=26, hjust = 0.025),
        axis.text=element_text(size=22),
        axis.text.x=element_text(angle=45,vjust=0.6),
        axis.title=element_text(size=26)
  )
```

## Tables:

```{r TableCorr}
CorrectTable$Variable<-as.character(CorrectTable$Variable)
CorrectTable$Variable<-gsub("MD","Ripple height",CorrectTable$Variable)
CorrectTable$Variable[CorrectTable$Variable=="scale(Start)"]<-"Day of experiment"
theme_zebra(regulartable(CorrectTable))%>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableCorr"), "Model output for discrimination performance"))

```
<br>
 
```{r TableSigmoid}
Sigtable$`Upper 95% CI`<-ifelse(Sigtable$`Upper 95% CI`>1000,9999,Sigtable$`Upper 95% CI`)
Sigtable$`Upper 95% CI`<-as.character(Sigtable$`Upper 95% CI`)
Sigtable$`Upper 95% CI`<-ifelse(Sigtable$`Upper 95% CI`=="9999",">1000",Sigtable$`Upper 95% CI`)
theme_zebra(regulartable(Sigtable)) %>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableSigmoid"), "Sigmoid fits of discrimination performance"))
```
<br>
```{r TableTime}
TimeTable$Variable<-as.character(TimeTable$Variable)
TimeTable$Variable<-gsub("MD","Ripple height",TimeTable$Variable)
TimeTable$Variable[TimeTable$Variable=="scale(Start)"]<-"Day of experiment"
theme_zebra(regulartable(TimeTable))%>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableTime"), "Trial time model outputs"))
```
<br>
```{r TableAbort}
AbortTable$Variable<-as.character(AbortTable$Variable)
AbortTable$Variable<-gsub("MD","Ripple height",AbortTable$Variable)
AbortTable$Variable[AbortTable$Variable=="scale(Start)"]<-"Day of experiment"
theme_zebra(regulartable(AbortTable)) %>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableAbort"), "Aborted trial model outputs"))
```
<br>
```{r TableDur}
DurTable$Variable<-as.character(DurTable$Variable)
DurTable$Variable<-gsub("mod_depth","Ripple height",DurTable$Variable)
DurTable$Variable[DurTable$Variable=="scale(start_t)"]<-"Day of experiment"
theme_zebra(regulartable(DurTable)) %>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableDur"), " sonar call duration model outputs"))

```
<br>
```{r TableRMS}
RMSTable$Variable<-as.character(RMSTable$Variable)
RMSTable$Variable<-gsub("mod_depth","Ripple height",RMSTable$Variable)
RMSTable$Variable[RMSTable$Variable=="start_t"]<-"Day of experiment"
theme_zebra(regulartable(RMSTable)) %>% 
  align(align="center",part="all") %>% 
  autofit%>% 
  add_header_lines(paste("Table", tab("fig: TableRMS"), " sonar call RMS model outputs"))

```
<br>
```{r TablePeak}
PeakTable$Variable<-as.character(PeakTable$Variable)
PeakTable$Variable<-gsub("mod_depth","Ripple height",PeakTable$Variable)
PeakTable$Variable[PeakTable$Variable=="scale(start_t)"]<-"Day of experiment"
theme_zebra(regulartable(PeakTable))%>% 
  align(align="center",part="all") %>% 
  autofit %>% 
  add_header_lines(paste("Table", tab("fig: TablePeak"), " peak frequency model outputs"))
```
<br>
```{r TableCent}
CentTable$Variable<-as.character(CentTable$Variable)
CentTable$Variable<-gsub("mod_depth","Ripple height",CentTable$Variable)
CentTable$Variable[CentTable$Variable=="scale(start_t)"]<-"Day of experiment"
theme_zebra(regulartable(CentTable))%>% 
  align(align="center",part="all") %>% 
  autofit %>% 
  add_header_lines(paste("Table", tab("fig: TableCent"), "Frequency centroid model outputs"))
```