##' Data Analysis of Electroantennogram Data from two different Stonefly Ecotypes
##' (full-winged and wing-reduced)
##'
##' Main R routine to reproduce the Fig 1, Fig 2, and Fig S1 from the manuscript
##' entitled 'Rapid evolution of olfactory degradation in recently flightless
##' insects'.
##' 
##

# load R packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(tibble)
library(ggbeeswarm)
library(scales)
library(arm)
library(blmeco)
library(lemon)
library(ggpubr)
library(grid)
library(gridExtra)
library(magick)
library(cowplot)
library(grImport)
library(pdftools)
library(multitaper)

# clean up environment
rm(list=ls()) 


# set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# load function to calculate various response dynamic measures
source("./functions/responseDynamics.R")
source('./functions/utils.R')

fig_path <- "../Figures"


### colors for stoneflies ###
fw_l <- "#787878" ## full-winged: grey
wr_l <- "#00CC00" ## wing-reduced: green


## load the sub plots for Fig. 1
Fig1A <- magick::image_read_pdf('../Pictures/Fig1A.pdf')
Fig1A <- ggplot() +
  theme_bw() +
  theme(panel.border = element_blank()) +
  draw_image(Fig1A)


Fig1B <- magick::image_read_pdf('../Pictures/Fig1B.pdf')
Fig1B <- ggplot() +
  theme_bw() +
  theme(panel.border = element_blank()) +
  draw_image(Fig1B)


# number of simulation draws from posterior distribution
nsim <- 10000

# aspect ratio (as y/x) for trace plots
ar <- 2/3

## Read data set
data2 <- read.csv('../Datasets/Dataset1.csv', stringsAsFactors = FALSE)


# create median filtered mean traces (only for non-10Hz stimulations)
med <- as.data.frame(t(apply(X = data2[data2$PulseDuration != 50, paste0("t", 1:5000)],
                            MARGIN = 1,
                            FUN = runmed, k = 11)))
data2[data2$PulseDuration != 50, paste0("t", 1:5000)] <- med


### automatically exclude alive Antennae that are not responding to 2-HeptanonTest ######
heptTest <- data2 %>%
  filter(Odor=="2HeptanoneTest", AntennaState=="alive")

par(mfrow=c(4,4))

analyze <- vector("logical", 0)
heptTest$keep <- NA

for(i in 1:nrow(heptTest)) {
  # in which ms is the signal higher than 2 x sd during 2 to 5 sec win of trace
  num <- which(heptTest[i, paste0("t", 250:1200)] > 2 * sd(heptTest[i, paste0("t", 2000:5000)]))
  result <- rle(diff(num))
  ## above threshold for at least 50 ms in a row?
  an <- any(result$lengths >= 50 & result$values==1)
  heptTest[i, "keep"] <- an 
  analyze <- c(analyze, an)
}

ID_keep <- unique(heptTest$ID_new[heptTest$keep])
ID_keep <- c(ID_keep, unique(data2[data2$Group=="PID", "ID_new"]))


data2 <- data2[data2$ID_new %in% ID_keep, ]


data2 <- droplevels(data2)


## Correct for amplification
## Note: 1 unit for different groups has different voltage ranges
## PID - 30 V
## Bee - 10 mV
## Fenestrata - 1 mV

multiply_voltage <- 1/data2$Amplification

data2[, paste0("t", 1:5000)] <- data2[, paste0("t", 1:5000)] * multiply_voltage # in Volt
data2[data2$Group != "PID", paste0("t", 1:5000)] <- data2[data2$Group != "PID", paste0("t", 1:5000)] * 1000 # in millivolt


## Convert some columns to a factor
data2$ID_new <- as.factor(data2$ID_new)
data2$Winged <- factor(data2$Winged, levels = c("full-winged", "wing-reduced", "none"))


data2$Location <- as.character(data2$Location)
data2$Location[data2$Group=="PID"] <- "PID"
data2$Location[data2$Species=="Honeybee"] <- "Bee"


# reorder the factor levels of column Odor and adjust the names
data2$Odor <- factor(data2$Odor, levels = c("Blank", "2Heptanone", "1Octanol", 
                                            "PropionicAcid", "2Butanone", "Stonefly", "2HeptanoneTest"),
                     labels = c("Blank","2-heptanone", "1-octanol", 
                                "Propionic acid", "2-butanone", "Stonefly", "2-heptanone (2nd)"))

data2$Group <- factor(data2$Group, levels = c("Burns",
                                              "BlackJack", "Whisky", "SixMile",
                                              "Lug", "Bee", "PID"),
                      labels = c("Mt Burns", "Black Jacks", "Whiskey", "Six Mile", "Lug", "Bee", "PID"))

data2 <- data2 %>%
  filter(AntennaState == "alive" | Group == "PID")


data2$NewGroup <- paste(data2$Group, data2$Winged)
data2$NewGroup <- as.factor(data2$NewGroup)



data2 <- droplevels(data2)


new_colnames <- c("Timeto10Max",
                  "Timeto90Max",
                  "TimetoMax", 
                  "TimeFromMaxto90",
                  "TimeFromMaxto10",
                  "Time10to90Max",
                  "Time90to10Max",
                  "Timeto10Min",
                  "Timeto90Min",
                  "TimetoMin", 
                  "TimeFromMinto90",
                  "TimeFromMinto10",
                  "Time10to90Min",
                  "Time90to10Min",
                  "MeanMax",
                  "MeanMin")

y1 <- as.data.frame(t(apply(data2[,paste0("t", 201:5000)], 1, responseDynamics, from=0.1, to=0.9)))

colnames(y1) <- new_colnames


data2 <- cbind(as.data.frame(data2), y1)

data2 <- data2 %>%
  mutate(Timeto10 = ifelse(Odor == "Propionic acid" |
                             Odor == "Stonefly", Timeto10Min, Timeto10Max),
         Time10to90 = ifelse(Odor == "Propionic acid" |
                               Odor == "Stonefly", Time10to90Min, Time10to90Max),
         Time90to10 = ifelse(Odor == "Propionic acid" |
                               Odor == "Stonefly", Time10to90Min, Time90to10Max))


data2$max <- apply(data2[,paste0("t", 251:1200)], 1, max)
data2$min <- apply(data2[,paste0("t", 251:1200)], 1, min)
data2 <- data2 %>%
  mutate(mima = ifelse(Odor == "Propionic acid" |
                         Odor == "Stonefly", min, max))
data2 <- data2 %>%
  mutate(Mean = ifelse(Odor == "Propionic acid" |
                         Odor == "Stonefly", MeanMin, MeanMax))

data2 <- data2 %>%
  mutate(TimeBackTo10 = ifelse(Odor == "Propionic acid" |
                                 Odor == "Stonefly", TimeFromMinto10, TimeFromMaxto10))




data2 %>%
  filter(Odor=="2-heptanone" | Odor=="2-heptanone (2nd)", PulseDuration==300) %>%
  group_by(Winged, ID_new, AntennaState) %>%
  summarise(ratio = mima[Odor=="2-heptanone (2nd)"] / mima[Odor=="2-heptanone"]) %>%
  group_by(Winged, AntennaState) %>%
  summarise(mean=mean(ratio), sd=sd(ratio), median=median(ratio), n=n())



#######################################
### How big is the sample size for the different groups?
data2 %>%
  group_by(Location, Group, Winged, ID_new) %>%
  summarise(n=n()) %>%
  group_by(Location, Group, Winged) %>%
  summarise(n=n())


tst <- data2
tst <- droplevels(tst)
table(tst$Config, tst$Group, tst$Winged)

locations <- c("Mt Burns", "Black Jacks", "Whiskey", "Six Mile", "Lug")
odors <- c("2-heptanone", "1-octanol", "2-butanone", "2-heptanone (2nd)")



############## Statistics #############
## Bayesian Statistics ###

creeks <- locations
plot_bayes_mima <- function(data, groups) {
  
  ################################################################################
  
  data1 <- data2
  
  odors <- c("2-heptanone", "1-octanol", "2-butanone")
  
  df <- as.data.frame(matrix(NA, ncol=7, nrow = 0))
  newDf <- as.data.frame(matrix(NA, ncol=3, nrow = 0))
  newDf_all <- as.data.frame(matrix(NA, ncol=3, nrow = 0))
  
  for (odor in odors) {
    if (odor == '2-heptanone') {
      data <- data1 %>% filter(Odor %in% c('2-heptanone (2nd)', '2-heptanone'))
    } else {
      data <- data1 %>% filter(Odor == odor)
    }
    
    data <- data %>% filter(Group %in% groups, PulseDuration != 50)
    
    data$PulseDuration <- as.character(data$PulseDuration)
    
    data$PulseDuration <- ifelse(data$PulseDuration == 300 & data$Odor == '2-heptanone (2nd)',
                                 '300\n (2nd)', data$PulseDuration)
    data$Odor <- as.character(data$Odor)
    
    data$Odor <- ifelse(data$Odor == '2-heptanone (2nd)', '2-heptanone', data$Odor)
    
    data$Location <- factor(data$Location)
    data$PulseDuration <- factor(data$PulseDuration, levels = c('15', '30', '150', '300', '300\n (2nd)'))
    data <- droplevels(data)
    
    
    
    
    for (cc in seq_along(groups)) {
      subdata <- data %>% filter(Group == groups[cc])
      subdata$PulseDuration <- as.factor(subdata$PulseDuration)
      subdata <- droplevels(subdata)
      min_mean <- min(subdata$Mean)
      
      if (min_mean < 0) {
        subdata$adj_Mean <- subdata$Mean + abs(min_mean) + 0.01
      } else {
        subdata$adj_Mean <- subdata$Mean
      }
      
      subdata$log2 <- log2(subdata$adj_Mean)
      
      if (nrow(subdata) == 0)
        next
      
      if (length(unique(subdata$PulseDuration)) == 1) {
        
        m1 <- lm(log2 ~ Winged, data=subdata)
        
      } else {
        
        if (unique(subdata$Winged == 'none')) {
          
          m1 <- lmer(log2 ~ PulseDuration + (1 | ID_new), data=subdata, REML = FALSE)
          
        } else {
          
          m1 <- lmer(log2 ~ Winged * PulseDuration + (1 | ID_new), data=subdata, REML = FALSE)
          
        }
        
        
      }
      
      
      par(mfrow=c(2,2))
      
      if (length(unique(subdata$PulseDuration)) != 1) {
        scatter.smooth(fitted(m1), resid(m1))
        abline(h=0, lty=2)
        title("Tukey-Anscombe Plot") # residuals vs. fitted
        
        qqnorm(resid(m1), main="normal QQ-plot, residuals") # qq of residuals
        qqline(resid(m1))
        scatter.smooth(fitted(m1), sqrt(abs(resid(m1)))) # res. var vs. fitted
        
        qqnorm(ranef(m1)$ID_new[,1], main="normal QQ-plot, random effects")
        qqline(ranef(m1)$ID_new[,1]) # qq of random effects
      } else {
        plot(m1)
      }
      
      
      summary(m1)      
      
      
      
      
      
      if (length(unique(subdata$PulseDuration)) == 1) {
        
        tmp <- unique(paste( subdata$PulseDuration, subdata$Winged, subdata$Odor, sep = '_'))
        
        tmp <- as.data.frame(matrix(unlist(strsplit(tmp, split = '_')), ncol=3, byrow = TRUE))
        
        tmp <- tmp[order(as.numeric(as.character(tmp$V1)), tmp$V2, tmp$V3), ]
        colnames(tmp) <- c('PulseDuration', 'Winged', 'Odor')
        
        tmp$PulseDuration <- factor(tmp$PulseDuration)
        tmp$PulseDuration <- factor(tmp$PulseDuration, 
                                    levels = levels(tmp$PulseDuration)[order(as.numeric(as.character(levels(tmp$PulseDuration))))])
        
        tmp$Winged <- factor(tmp$Winged, levels = c('full-winged', 'wing-reduced'))
        tmp$Odor <- factor(tmp$Odor, levels = c('2-heptanone', '1-octanol', '2-butanone'))
        
        newdat <- tmp
        
        bsim <- sim(m1, n.sim=nsim)
        
        fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
        Xmat <- model.matrix( ~ Winged, data=newdat)
        
        
        for(i in 1:nsim) fitmat[,i] <- Xmat %*% bsim@coef[i,] #fitted values
        
        
        tmp <- t(apply(fitmat, 1, quantile, prob=c(0.025, 0.5, 0.975)))
        
        if (min_mean < 0) {
          
          tmp <- 2^(tmp) - abs(min_mean) - 0.01
          
        } else {
          
          tmp <- 2^(tmp)
          
        }
        
        
        
        
        tst <- cbind(newdat, tmp, 'Group' = groups[cc])
        
        df <- rbind(df, tst)
        
        for (oo in levels(newdat$Odor)) {
          
          datt <- newdat[newdat$Odor == oo, ]
          
          for (pp in unique(datt$PulseDuration)) {
            
            dat <- fitmat[which(newdat$PulseDuration == pp & newdat$Odor == oo), ]
            star2 <- sum(dat[1, ] > dat[2, ]) / nsim  
            
            newDf <- rbind(newDf, 
                           data.frame('Group' = groups[cc],
                                      'Odor' = oo,
                                      'Diff' = star2,
                                      'PulseDuration' = pp))
            
          }
          
          
        }
        
        
      } else {
        
        if (unique(subdata$Winged == 'none')) {
          
          #######################
          
          tmp <- unique(paste(subdata$PulseDuration, subdata$Winged, subdata$Odor, sep = '_'))
          
          tmp <- as.data.frame(matrix(unlist(strsplit(tmp, split = '_')), ncol=3, byrow = TRUE))
          
          tmp <- tmp[order(as.numeric(as.character(tmp$V1)), tmp$V2, tmp$V3), ]
          colnames(tmp) <- c('PulseDuration', 'Winged', 'Odor')
          
          tmp$PulseDuration <- factor(tmp$PulseDuration)
          tmp$PulseDuration <- factor(tmp$PulseDuration, 
                                      levels = levels(tmp$PulseDuration)[order(as.numeric(as.character(levels(tmp$PulseDuration))))])
          tmp$Winged <- factor(tmp$Winged, levels = c('full-winged', 'wing-reduced', 'none'))
          tmp$Odor <- factor(tmp$Odor, levels = c('2-heptanone', '1-octanol', '2-butanone'))
          
 
          newdat <- tmp
          
          
          bsim <- sim(m1, n.sim=nsim)
          
          fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
          Xmat <- model.matrix( ~ PulseDuration, data=newdat)
          
          
          for(i in 1:nsim) fitmat[,i] <- Xmat %*% bsim@fixef[i,] #fitted values
          
          names(Xmat)
          length(newdat$Odor)
          
          tmp <- t(apply(fitmat, 1, quantile, prob=c(0.025, 0.5, 0.975)))
          
          if (min_mean < 0) {
            
            tmp <- 2^(tmp) - abs(min_mean) - 0.01
            
          } else {
            
            tmp <- 2^(tmp)
            
          }
          
          tst <- cbind(newdat, tmp, 'Group' = groups[cc])
          
          df <- rbind(df, tst)
          
          
          for (oo in levels(newdat$Odor)) {
            
            datt <- newdat[newdat$Odor == oo, ]
            
            for (pp in unique(datt$PulseDuration)) {
              
              dat <- fitmat[which(newdat$PulseDuration == pp & newdat$Odor == oo), ]

              newDf <- rbind(newDf, 
                             data.frame('Group' = groups[cc],
                                        'Odor' = oo,
                                        'Diff' = 0.5,
                                        'PulseDuration' = pp))
              
            }
            
            
          }
          
          #######################
          
        } else {
          
          tmp <- unique(paste( subdata$PulseDuration, subdata$Winged, subdata$Odor, sep = '_'))
          
          tmp <- as.data.frame(matrix(unlist(strsplit(tmp, split = '_')), ncol=3, byrow = TRUE))
          
          tmp <- tmp[order(as.numeric(as.character(tmp$V1)), tmp$V2, tmp$V3), ]
          colnames(tmp) <- c('PulseDuration', 'Winged', 'Odor')
          
          tmp$PulseDuration <- factor(tmp$PulseDuration)
          tmp$PulseDuration <- factor(tmp$PulseDuration, 
                                      levels = levels(tmp$PulseDuration)[order(as.numeric(as.character(levels(tmp$PulseDuration))))])
          tmp$Winged <- factor(tmp$Winged, levels = c('full-winged', 'wing-reduced'))
          tmp$Odor <- factor(tmp$Odor, levels = c('2-heptanone', '1-octanol', '2-butanone'))
          
          
          newdat <- tmp
          
          
          bsim <- sim(m1, n.sim=nsim)
          
          fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
          Xmat <- model.matrix( ~ Winged * PulseDuration, data=newdat)
          
          
          
          for(i in 1:nsim) fitmat[,i] <- Xmat %*% bsim@fixef[i,] #fitted values
          
          tmp <- t(apply(fitmat, 1, quantile, prob=c(0.025, 0.5, 0.975)))
          
          if (min_mean < 0) {
            
            tmp <- 2^(tmp) - abs(min_mean) - 0.01
            
          } else {
            
            tmp <- 2^(tmp)
            
          }
          
          tst <- cbind(newdat, tmp, 'Group' = groups[cc])
          
          df <- rbind(df, tst)
          
          
          for (oo in levels(newdat$Odor)) {
            
            datt <- newdat[newdat$Odor == oo, ]
            
            for (pp in unique(datt$PulseDuration)) {
              
              dat <- fitmat[which(newdat$PulseDuration == pp & newdat$Odor == oo), ]
              star2 <- sum(dat[1, ] > dat[2, ]) / nsim  
              
              newDf <- rbind(newDf, 
                             data.frame('Group' = groups[cc],
                                        'Odor' = oo,
                                        'Diff' = star2,
                                        'PulseDuration' = pp))
              
            }
            
          }
          
        }
        
      }
      
    }
    
  }
  
  
  colnames(df)[4:6] <- c("lower", "fitted", "upper")
  newDf
  
  newDf$Star1 <- newDf$Diff >= 0.95 | newDf$Diff <= 0.05
  newDf$Star2 <- newDf$Diff >= 0.99 | newDf$Diff <= 0.01
  newDf$Star <- ''
  newDf$Star[newDf$Star1 == TRUE] <- '*'
  newDf$Star[newDf$Star2 == TRUE] <- '**'
  
  df$Group
  newDf$Winged <- 'star'

  col <- c(fw_l, wr_l, "blue")
  
  plotdata <- data1 %>%
    filter(Odor %in% c("2-heptanone", "2-heptanone (2nd)"), 
           PulseDuration != 50)
  plotdata$PulseDuration <- as.character(plotdata$PulseDuration)
  plotdata$PulseDuration[plotdata$Odor == "2-heptanone (2nd)"] <- '300\n (2nd)'
  
  plotdata$Odor[plotdata$Odor == "2-heptanone (2nd)"] <- "2-heptanone"
  
  plotdata$PulseDuration <- factor(plotdata$PulseDuration, levels = c(15, 30, 150, 300, "300\n (2nd)"))
  
  
  plotdf <- df %>%
    filter(Odor %in% c("2-heptanone", "2-heptanone (2nd)"))
  plotdf$PulseDuration <- as.character(plotdf$PulseDuration)
  plotdf$PulseDuration[plotdf$Odor == "2-heptanone (2nd)"] <- '300\n (2nd)'
  plotdf$Odor[plotdf$Odor == "2-heptanone (2nd)"] <- "2-heptanone"
  plotdf$PulseDuration <- factor(plotdf$PulseDuration, levels = c(15, 30, 150, 300, "300\n (2nd)"))
  
  
  newDf_plot <- newDf %>%
    filter(Odor %in% c("2-heptanone", "2-heptanone (2nd)"))
  newDf_plot$PulseDuration <- as.character(newDf_plot$PulseDuration)
  newDf_plot$PulseDuration[newDf_plot$Odor == "2-heptanone (2nd)"] <- '300\n (2nd)'
  newDf_plot$Odor[newDf_plot$Odor == "2-heptanone (2nd)"] <- "2-heptanone"
  newDf_plot$PulseDuration <- factor(newDf_plot$PulseDuration, levels = c(15, 30, 150, 300, "300\n (2nd)"))
  
  
  fen <- ggplot(plotdata %>% filter(Group %in% creeks), aes(PulseDuration, Mean, color=Winged)) +
    facet_grid(Odor ~ Group) +
    geom_beeswarm(dodge.width = 0.5, pch = 1) +
    geom_errorbar(data=plotdf %>% filter(Group %in% creeks), aes(ymin=lower, ymax=upper, y=NULL),
                  color="black", width=0.5, position = position_dodge2(padding = 1, width=0.1)) +
    stat_summary(data=newDf_plot %>% filter(Group %in% creeks),
                 aes(x=as.factor(PulseDuration), y=0.15, label = Star),
                 geom = 'text', fun = max, color="black",
                 size = 8) +
    geom_beeswarm(data=plotdf %>% filter(Group %in% creeks), aes(y = fitted, group = Winged),
                  pch = '_', dodge.width = 0.5, color = 'black', size=3) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    theme_bw() +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(colour = "black"),
          axis.title = element_blank(),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          text = element_text(size=16),
          axis.text = element_text(colour="black"))
  fen
  
  bee <- ggplot(plotdata %>% filter(Group %in% 'Bee'), aes(PulseDuration, Mean, color=Winged)) +
    facet_grid(Odor ~ Group) +
    geom_beeswarm(dodge.width = 0.5, pch = 1) +
    geom_errorbar(data=plotdf %>% filter(Group %in% 'Bee'), aes(ymin=lower, ymax=upper, y=NULL), color="black",
                  width=0, position = position_dodge2(padding = 1, width=0.1)) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = 'blue') +
    theme_bw() +
    geom_beeswarm(data=plotdf %>% filter(Group %in% 'Bee'), aes(y=fitted, group=Winged),
                  pch = '_', color = 'black', groupOnX = TRUE, size=3) +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(colour = "black"),
          axis.title=element_blank(),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          text = element_text(size=16),
          axis.text=element_text(colour="black"))
  bee
  
  pid <- ggplot(plotdata %>% filter(Group %in% 'PID'), aes(PulseDuration, Mean, color=Winged)) +
    facet_grid(Odor ~ Group) +
    geom_beeswarm(dodge.width = 0.5, pch = 1) +
    geom_errorbar(data=plotdf %>% filter(Group %in% 'PID'), aes(ymin=lower, ymax=upper, y=NULL),
                  color="black",  width=0, position = position_dodge2(padding = 1, width=0.1)) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = 'blue') +
    theme_bw() +
    geom_beeswarm(data=plotdf %>% filter(Group %in% 'PID'), aes(y=fitted, group=Winged),
                  pch = '_', dodge.width = 0.5, color = 'black',  groupOnX = TRUE, size=3) +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(colour = "black"),
          axis.title=element_blank(),
          plot.tag = element_text(size = 14, hjust = 0.5, vjust = 0.5),
          plot.tag.position = c(0.02, 1 - 0.03),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          text = element_text(size=16),
          axis.text=element_text(colour="black"))
  pid
  
  l <- list('Diffs' = newDf, 'Mean_crI' = df, 'fen' = fen, 'bee' = bee, 'pid' = pid)
  return(l)
}


### Figure S1 #####
Fly_mean <- plot_bayes_mima(data = data2, groups = c(creeks, 'Bee', 'PID'))

FigS1B <- egg::ggarrange(Fly_mean$fen, Fly_mean$bee, Fly_mean$pid, clip = FALSE,
                      labels = c('', '', expression(x ~ 10^3)), label.args = list(vjust=0),
                      nrow = 1, ncol=3, widths = c(5, 1, 1), heights = c(6), debug = FALSE)

y.grob <- textGrob(paste("Response strength (mV)"), 
                   gp=gpar(fontsize=18), rot=90)

x.grob <- textGrob("Pulse duration (ms)", 
                   gp=gpar(fontsize=18))

#add to plot
FigS1B <- grid.arrange(arrangeGrob(FigS1B, left = y.grob, bottom = x.grob))



plot_winged_all <- function(data, odors, creeks, Diffs) {
  
  data <- data2
  data <- data %>%
    filter(Group %in% c(creeks, "Bee", "PID")) %>%
    filter(PulseDuration == 300) %>%
    filter(Odor %in% odors)
  
  data <- droplevels(data)
  
  data$Odor <- factor(data$Odor, levels = levels(data$Odor),
                      labels = c("2-heptanone", "1-octanol", "2-butanone", "2-heptanone\n (2nd)"))
  
  # calculate the mean over all odor-pulseduration configurations
  d_mean <- data %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise_at(.funs = mean, .vars=paste0("t", 1:5000))
  
  
  # calculate the standard error over all odor-pulseduration configurations
  d_se <- data %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise_at(.funs = se, .vars=paste0("t", 1:5000))
  
  # calculate sample size (n) over all odor-pulseduration configurations
  d_n <- data %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise(N=n()) %>%
    group_by(Group, Winged) %>%
    summarise(N=unique(N))
  
  d_max <- d_mean %>%
    group_by(Group, Odor, PulseDuration) %>%
    summarise(Max=max(t400))
  
  
  # create d_all and add the sample size
  d_all <- d_mean
  data_se <- gather(d_se, time, se, t1:t5000, factor_key=TRUE)
  
  
  d_all$PulseDuration <- factor(d_all$PulseDuration, levels = c(300))
  d_all$PulseDur <-  factor(d_all$PulseDuration, levels = c(300),
                            labels = c("300 ms"))
  
  
  data_long <- gather(d_all, time, Voltage, t1:t5000, factor_key=TRUE)
  data_long$time <- (as.numeric(gsub("t", "", data_long$time)) - 200)/1000
  data_long$SE <- data_se$se
  
  
  data_long$high <- data_long$Voltage + data_long$SE
  data_long$low <- data_long$Voltage - data_long$SE
  
  
  tmp <- expand.grid(creeks, c("FW", "WR"), stringsAsFactors = FALSE)
  tmp <- tmp[order(tmp$Var1), ]
  tmp <- paste0(tmp$Var2, ", ", tmp$Var1) 
  
  tmp <- c(tmp, "Bee", "PID")
  
  data_long$NewGroup <- factor(paste(data_long$Winged, data_long$Group),
                               labels = tmp)
  
  col <- c(fw_l, wr_l, "blue")
  
  
  print(d_n)
  
  geom.text.size <- 5.2
  data_long$Winged <- factor(data_long$Winged,
                             levels = levels(data_long$Winged),
                             labels = c("Full-winged", "Wing-reduced", "none"))
  
  
  data_segm1 <- data.frame(y=0, yend=0.1, x=Inf, xend=Inf,
                           PulseDuration= c(300),
                           Group="Lug",
                           Odor=c("2-heptanone\n (2nd)"))
  label1 <- data.frame(x=Inf, y=0.05,
                       text=c("0.1 mV"),
                       PulseDuration= c(300),
                       Group="Lug",
                       Odor=c("2-heptanone\n (2nd)"))
  
  fw_text <- data.frame(x = 1, y = 0.07,
                        text=c("full-winged"),
                        PulseDuration= c(300),
                        Group="Mt Burns",
                        Odor=c("2-heptanone"))
  wr_text <- data.frame(x = 1, y = 0.1,
                        text=c("wing-reduced"),
                        PulseDuration= c(300),
                        Group="Mt Burns",
                        Odor=c("2-heptanone"))
  
  Diffs$Odor <- as.character(Diffs$Odor)
  Diffs$Odor[Diffs$PulseDuration == '300\n (2nd)'] <- '2-heptanone\n (2nd)'
  Diffs$Odor <- factor(Diffs$Odor)
  Diffs$Odor <- factor(Diffs$Odor, levels = c('2-heptanone', '1-octanol', '2-butanone', '2-heptanone\n (2nd)'))
  
  Diffs$PulseDuration[Diffs$PulseDuration == '300\n (2nd)'] <- '300'
  
  Diffs <- droplevels(Diffs)
  
  Diffs <- merge(Diffs, d_max, by = c('Group', 'Odor', 'PulseDuration'))
  
  ## plotting mean and se
  Fen <- ggplot(data_long %>% filter(Group %in% creeks), aes(x=time, y=Voltage)) + 
    facet_grid(Odor ~ Group, switch = "y") +
    geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
    geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    theme_bw() +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_blank(),
          axis.title = element_blank(),
          plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          strip.background = element_rect(colour="NA", fill="NA"),
          panel.border = element_blank(),
          text = element_text(size=18),
          axis.ticks = element_blank(),
          aspect.ratio = ar * 4,
          strip.text.y.left = element_text(angle = 0, size = 16),
          strip.text.x = element_text(size = 16),
          axis.text=element_blank()) +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
    geom_text(data = wr_text, aes(x=x, y=y, label=text), size=4, hjust=0, color=wr_l) +
    geom_text(data = fw_text, aes(x=x, y=y, label=text), size=4, hjust=0, color=fw_l) +
    geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
    coord_cartesian(expand = TRUE,
                    clip = 'off')
  Fen
  
  data_segm1 <- data.frame(y=0, yend=1, x=Inf, xend=Inf,
                           PulseDuration= c(300),
                           Group="Bee",
                           Odor=c("2-heptanone\n (2nd)"))
  label1 <- data.frame(x=Inf, y=0.5,
                       text=c("1 mV"),
                       PulseDuration= c(300),
                       Group="Bee",
                       Odor=c("2-heptanone\n (2nd)"))
  seg <- data.frame(x=4.5, y = -0.35,
                    text=c(""),
                    PulseDuration= c(300),
                    Group="Bee",
                    Odor=c("2-heptanone\n (2nd)"))
  
  
  B <- ggplot(data_long %>% filter(Group %in% "Bee"), aes(x=time, y=Voltage)) + 
    facet_grid(Odor ~ Group, switch = "y") +
    geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = "Blue") +
    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
    geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    theme_bw() +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_blank(),
          axis.title = element_blank(),
          plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          strip.background = element_rect(colour="NA", fill="NA"),
          panel.border = element_blank(),
          text = element_text(size=18),
          axis.ticks = element_blank(),
          aspect.ratio = ar * 4,
          strip.text.y.left = element_blank(),
          strip.text.x = element_text(size = 16),
          axis.text=element_blank()) +
    geom_blank(data = seg, mapping = aes(x=x, y=y), inherit.aes = FALSE) +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
    coord_cartesian(expand = TRUE,
                    clip = 'off')
  
  data_segm1 <- data.frame(y=0, yend=c(3, 0.05, 3), x=Inf, xend=Inf,
                           PulseDuration= c(300),
                           Group="PID",
                           Odor=c("2-heptanone\n (2nd)", "1-octanol", "2-heptanone"))
  label1 <- data.frame(x=Inf, y=c(1.5, 0.04, 1.5),
                       text=c("3 V", "0.05 V", "3 V"),
                       PulseDuration= c(300),
                       Group="PID",
                       Odor=c("2-heptanone\n (2nd)",  "1-octanol","2-heptanone"))
  
  data_segm2 <- data.frame(x=3.5, xend=4.5, y=-0.3, yend=-0.3,
                           PulseDuration= c(300),
                           Group="PID",
                           Odor=c("2-heptanone\n (2nd)"))
  
  label2 <- data.frame(x=4, y=-0.2,
                       text=c("1 s"),
                       PulseDuration= c(300),
                       Group="PID",
                       Odor=c("2-heptanone\n (2nd)"))
  
  seg <- data.frame(x=4.5, y = -0.35,
                    text=c(""),
                    PulseDuration= c(300),
                    Group="PID",
                    Odor=c("2-butanone"))
  
  P <- ggplot(data_long %>% filter(Group %in% "PID"), aes(x=time, y=Voltage)) + 
    facet_grid(Odor ~ Group, switch = "y", scales = 'free_y') +
    geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = "Blue") +
    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
    geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    theme_bw() +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_blank(),
          axis.title = element_blank(),
          plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          strip.background = element_rect(colour="NA", fill="NA"),
          panel.border = element_blank(),
          text = element_text(size=18),
          axis.ticks = element_blank(),
          aspect.ratio = ar * 4,
          strip.text.y.left = element_blank(),
          strip.text.x = element_text(size = 16),
          axis.text=element_blank()) +
    geom_blank(data = seg, mapping = aes(x=x, y=y), inherit.aes = FALSE) +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
    geom_segment(data=data_segm2,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label2, aes(x=x, y=y, label=text), size=geom.text.size, vjust=1.3) +
    coord_cartesian(expand = TRUE,
                    clip = 'off')
  P
  
  cowplot::plot_grid(Fen, B, P, nrow = 1, align = "h")
  list1 <- list("traceW"=Fen,
                "traceB"=B,
                "traceP"=P)
}


traces300 <- plot_winged_all(data = data2,
                             odors = odors,
                             creeks = locations,
                             Diffs = Fly_mean$Diffs %>% filter(Group %in% creeks,
                                                               PulseDuration %in% c('300', '300\n (2nd)')))

Fig1C <- egg::ggarrange(traces300$traceW, traces300$traceB, traces300$traceP,
                        heights = 5, nrow = 1, ncol=3, widths = c(5, 1, 1), debug = FALSE)


########################

data <- data2
data <- data %>%
  filter(Group %in% c(creeks, "Bee", "PID")) %>%
  filter(Odor %in% c("2-heptanone"), !Frequency %in% 10, !PulseDuration %in% 300)

data <- droplevels(data)

data$Odor <- factor(data$Odor, levels = levels(data$Odor),
                    labels = c("2-heptanone"))

# calculate the mean over all odor-pulseduration configurations
d_mean <- data %>%
  group_by(Group, Winged, Odor, PulseDuration) %>%
  summarise_at(.funs = median, .vars=paste0("t", 1:5000))


# calculate sample size (n) over all odor-pulseduration configurations
d_n <- data %>%
  group_by(Group, Winged, Odor, PulseDuration) %>%
  summarise(N=n()) %>%
  group_by(Group, Winged) %>%
  summarise(N=unique(N))

d_max <- d_mean %>%
  group_by(Group, Odor, PulseDuration) %>%
  summarise(Max=max(t400))


# create d_all and add the sample size
d_all <- d_mean

d_all$PulseDuration <- factor(d_all$PulseDuration, levels = c('15', '30', '150'))
d_all$PulseDur <-  factor(d_all$PulseDuration)


data_long <- gather(d_all, time, Voltage, t1:t5000, factor_key=TRUE)
data_long$time <- (as.numeric(gsub("t", "", data_long$time)) - 200)/1000





print(d_n)

geom.text.size <- 5.2
data_long$Winged <- factor(data_long$Winged,
                           levels = levels(data_long$Winged),
                           labels = c("Full-winged", "Wing-reduced", "none"))


data_segm1 <- data.frame(y=c(0, 0), yend=c(0.1, 0), x=c(Inf, 2), xend=c(Inf, 3),
                         PulseDuration = factor(150, levels = c(15, 30, 150)),
                         Group=factor(c("Lug", 'Mt Burns'), levels = creeks),
                         Odor=c("2-heptanone"))
label1 <- data.frame(x=c(Inf, 3), y=c(0.05, 0.02),
                     text=c("0.1 mV", '1 s'),
                     PulseDuration= factor(150, levels = c(15, 30, 150)),
                     Group=factor(c("Lug", 'Mt Burns'), levels = creeks),
                     Odor=c("2-heptanone"))

fw_text <- data.frame(x = 0, y = 0.02,
                      text=c("full-winged"),
                      PulseDuration= factor(15, levels = c(15, 30, 150)),
                      Group=factor("Mt Burns", levels = creeks),
                      Odor=c("2-heptanone"))
wr_text <- data.frame(x = 0, y = 0.05,
                      text=c("wing-reduced"),
                      PulseDuration = factor(15, levels = c(15, 30, 150)),
                      Group=factor("Mt Burns", levels = creeks),
                      Odor=c("2-heptanone"))

pulseDur_text <- data.frame(x = 0.25, y = 0.1,
                      text=c("15 ms", "30 ms", "150 ms"),
                      PulseDuration = factor(c("15", "30", "150"), levels = c("15", "30", "150")),
                      Group=factor("Black Jacks", levels = creeks),
                      Odor=c("2-heptanone"))

data_long$Group <- factor(data_long$Group, levels = c(creeks, 'Bee', 'PID'))


col <- c(fw_l, wr_l)

Fen <- ggplot(data_long %>% filter(Group %in% creeks), aes(x=time, y=Voltage)) + 
  facet_grid(PulseDuration ~ Group, switch = "y") +
  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                        ymin=-Inf, ymax=Inf),
            fill="grey90", alpha=1) +
  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
  geom_line(aes(colour = Winged)) +
  theme_bw() +
  labs(y = ' ') +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        panel.border = element_blank(),
        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
        strip.background = element_rect(colour="NA", fill="NA"),
        text = element_text(size=18),
        aspect.ratio = ar * 3,
        strip.text.y.left = element_blank(),
        strip.text.x = element_text(size = 20),
        axis.text=element_blank()) +
  geom_segment(data=data_segm1,
               aes(x=x, y=y, yend=yend, xend=xend)) +
  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
  geom_text(data = wr_text, aes(x=x, y=y, label=text), size=6, hjust=0, color=wr_l) +
  geom_text(data = fw_text, aes(x=x, y=y, label=text), size=6, hjust=0, color=fw_l) +
  geom_text(data = pulseDur_text, aes(x=x, y=y, label=text), size=5, hjust=0, color='black') +
  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
  coord_cartesian(expand = TRUE,
                  clip = 'off')
Fen

col <- c('blue')
data_segm1 <- data.frame(y=0, yend=1, x=Inf, xend=Inf,
                         PulseDuration = factor(150, levels = c(15, 30, 150)),
                         Group=factor("Bee"),
                         Odor=c("2-heptanone"))
label1 <- data.frame(x=Inf, y=0.5,
                     text=c("1 mV"),
                     PulseDuration= factor(150, levels = c(15, 30, 150)),
                     Group=factor("Bee"),
                     Odor=c("2-heptanone"))


Bee <- ggplot(data_long %>% filter(Group %in% 'Bee'), aes(x=time, y=Voltage)) + 
  facet_grid(PulseDuration ~ Group, switch = "y") +
  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                        ymin=-Inf, ymax=Inf),
            fill="grey90", alpha=1) +
  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
  geom_line(aes(colour = Winged)) +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
        strip.background = element_rect(colour="NA", fill="NA"),
        panel.border = element_blank(),
        text = element_text(size=18),
        axis.ticks = element_blank(),
        aspect.ratio = ar * 3,
        strip.text.y = element_blank(),
        strip.text.x = element_text(size = 20),
        axis.text=element_blank()) +
  geom_segment(data=data_segm1,
               aes(x=x, y=y, yend=yend, xend=xend)) +
  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
  coord_cartesian(expand = TRUE,
                  clip = 'off')
Bee

data_segm1 <- data.frame(y=0, yend=3, x=Inf, xend=Inf,
                         PulseDuration = factor(150, levels = c(15, 30, 150)),
                         Group=factor("PID"),
                         Odor=c("2-heptanone"))
label1 <- data.frame(x=Inf, y=1.5,
                     text=c("3 V"),
                     PulseDuration= factor(150, levels = c(15, 30, 150)),
                     Group=factor("PID"),
                     Odor=c("2-heptanone"))

PID <- ggplot(data_long %>% filter(Group %in% 'PID'), aes(x=time, y=Voltage)) + 
  facet_grid(PulseDuration ~ Group, switch = "y") +
  geom_rect(mapping=aes(xmin=0, xmax=as.numeric(as.character(PulseDuration))/1000,
                        ymin=-Inf, ymax=Inf),
            fill="grey90", alpha=1) +
  scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
  geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
  geom_line(aes(colour = Winged)) +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(5.5, 5.5, 5.5, 5.5), "points"),
        strip.background = element_rect(colour="NA", fill="NA"),
        panel.border = element_blank(),
        text = element_text(size=18),
        axis.ticks = element_blank(),
        aspect.ratio = ar * 3,
        strip.text.y = element_blank(),
        strip.text.x = element_text(size = 20),
        axis.text=element_blank()) +
  geom_segment(data=data_segm1,
               aes(x=x, y=y, yend=yend, xend=xend)) +
  geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
  # geom_text(data = Diffs, aes(x=0.4, y=Max + 0.025, label=Star), size=geom.text.size, hjust=0.5) +
  coord_cartesian(expand = TRUE,
                  clip = 'off')
PID


FigS1A <- egg::ggarrange(Fen, Bee, PID,
                         heights = 5, nrow = 1, ncol=3, widths = c(5, 1, 1), debug = FALSE)

########################


#### Function to plot parameters for Response Dynamic:

plot_bayes_time <- function(data, creeks, yvalue, ylab) {
  
  data <- data2
  
  data <- data %>%
    filter(PulseDuration==300, Group %in% creeks, Odor %in% odors)
  
  df <- as.data.frame(matrix(NA, ncol=7, nrow = 0))
  df_diff <- as.data.frame(matrix(NA, ncol=8, nrow = 0))
  
  ##########################################################################
  
  for (cc in seq_along(creeks)) {
    
    subdata <- data %>%
      filter(Group == creeks[cc])
    
    subdata$PulseDuration <- as.factor(subdata$PulseDuration)
    subdata <- droplevels(subdata)
    
    m1 <- lm(log2(get(yvalue)) ~ Winged * Odor, data = subdata)
    par(mfrow=c(2,2))
    plot(m1)
    summary(m1)
    
    subdata$Odor
    
    newdat <- expand.grid(Winged = factor(c("full-winged","wing-reduced")),
                          PulseDuration = factor(levels(subdata$PulseDuration), levels = levels(subdata$PulseDuration)),
                          Odor = factor(levels(subdata$Odor), levels = odors))
    
    bsim <- sim(m1, n.sim=nsim)
    
    fitmat <- matrix(ncol=nsim, nrow=nrow(newdat))
    Xmat <- model.matrix( ~ Winged * Odor, data=newdat)
    for(i in 1:nsim) fitmat[,i] <- Xmat %*% bsim@coef[i,] #fitted values
    
    tmp <- t(apply(fitmat, 1, quantile, prob=c(0.025, 0.5, 0.975)))
    
    tst <- cbind(newdat, tmp, 'Group' = creeks[cc])
    
    df <- rbind(df, tst)
    
    for (oo in odors) {
      
      dat <- fitmat[which(newdat$PulseDuration == 300 & newdat$Odor == oo), ]
      star2 <- sum(dat[1, ] > dat[2, ]) / nsim  
      
      
      
      df_diff <- rbind(df_diff, 
                       data.frame('Group' = creeks[cc],
                                  'Odor' = oo,
                                  'Diff' = star2))
      
    }
    
    
  }
  
  df_diff
  df_diff$Star1 <- df_diff$Diff >= 0.95 | df_diff$Diff <= 0.05
  df_diff$Star2 <- df_diff$Diff >= 0.99 | df_diff$Diff <= 0.01
  
  
  colnames(df)[4:6] <- c('lower', 'fitted', 'upper')
  
  
  df$lower <- 2^(df$lower)
  df$fitted <- 2^(df$fitted)
  df$upper <- 2^(df$upper)
  
  
  subdf <- df %>%
    group_by(Group, Odor) %>%
    summarise(max=max(fitted))
  
  
  sub_data2 <- data %>%
    filter(Group!="Fenestrata-Bee",
           Group!="Bee",
           Odor!="Blank",
           Odor!="Propionic acid",
           Odor!="Stonefly")
  
  sub_data2 <- droplevels(sub_data2)
  df <- droplevels(df)
  df_diff <- droplevels(df_diff)
  
  sub_data2$PulseDuration <- factor(sub_data2$PulseDuration,
                                    levels = c(300))
  
  df$PulseDuration <- factor(df$PulseDuration,
                             levels = c(300))
  
  
  df$NewGroup <- as.factor(paste(df$Winged, df$Group))
  
  df_diff$NewGroup <- as.factor(paste(df_diff$Winged, df_diff$Group))
  
  
  
  col <- c(fw_l, wr_l)
  col <- c(col[1], "black", col[2])
  
  df_diff$Star[df_diff$Star1 == TRUE] <- "*"
  df_diff$Star[df_diff$Star1 == FALSE] <- " "
  df_diff$Star[df_diff$Star2 == TRUE] <- "**"
  stars <- df_diff$Star
  
  df_diff$NewGroup <- "star"
  
  print(df_diff)
  
  df_diff$Winged <- "star"
  
  df$Odor <- factor(df$Odor, levels = odors)
  df_diff$Odor <- factor(df_diff$Odor, levels = odors)
  sub_data2$Odor <- factor(sub_data2$Odor, levels = odors)
  
  MAX <- sub_data2 %>%
    group_by(Odor, Group) %>%
    summarise(max = max(get(yvalue), na.rm = TRUE))
  df_diff <- merge(df_diff, MAX, by=c('Odor', 'Group'))
  
  p <- ggplot(df, aes(x=Group, y=fitted, color=Winged)) +
    geom_hline(aes(yintercept=0), color="black", size=0.3, linetype="dotted") +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    scale_y_continuous(breaks = pretty_breaks(n=2), expand = expansion(mult = c(0.05, .1))) +
    facet_grid(Odor ~ Group, scales = "free") +
    geom_beeswarm(sub_data2, dodge.width = 0.8, mapping=aes(x=Group, y=get(yvalue)), cex=1, shape=21) +
    geom_point(aes(group=NewGroup), position = position_dodge(0.8), shape="_", cex=4, color="black") +
    geom_errorbar(aes(ymin=lower, ymax=upper, group=Winged),
                  color="black", width=0, position=position_dodge(0.8), size=0.5) +
    theme_bw() +
    ylab(ylab) +
    xlab(" ") +
    scale_x_discrete(labels=c("", "")) +
    geom_text(data=df_diff, aes(x=Group, y = max * 0.9, label=Star), size = 7) +
    theme(legend.position = "none",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line.x = element_blank(),
          axis.line.y = element_line(colour = "black"),
          axis.title.x=element_blank(),
          axis.ticks.x=element_blank(),
          strip.text.x = element_text(size = 16),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text.y = element_text(size = 16, angle = 0, hjust = 0),
          panel.border = element_blank(),
          text = element_text(size=18),
          axis.text=element_text(colour="black"))
  print(p)
  
  return(p)
}



####### Figure 2 ##########################

Fig2A <- plot_bayes_time(data = data2, yvalue = "Timeto10", ylab = "Response onset (ms)", creeks = creeks)
Fig2B <- plot_bayes_time(data2, "TimeBackTo10", "Response offset (ms)", creeks = creeks)


##################  10 Hz Stimulations ###########
on_times <- seq(0, 2.9, by = 0.1)
off_times <- seq(0.05, 3, by=0.1)

times_tst_on <- expand.grid(levels(data2$Odor), "on"=on_times)
times_tst_off <- expand.grid(levels(data2$Odor), "off"=off_times)

Times <- data.frame("Odor"=times_tst_on$Var1, "on"=times_tst_on$on, "off"=times_tst_off$off)

Times <- Times[!(Times$Odor=="Blank"), ]

plot_freq <- function(data, group, times) {
  
  data <- data %>%
    filter(PulseDuration == 50,
           Odor == "2-heptanone")
  
  data <- droplevels(data)
  
  d_mean <- data %>%
    filter(Group %in% group) %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise_at(.funs = median, .vars=paste0("t", 1:5000))
  
  # calculate the standard error 
  d_se <- data %>%
    filter(Group %in% group) %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise_at(.funs = se, .vars=paste0("t", 1:5000))
  
  # calculate sample size (n) 
  d_n <- data %>%
    filter(Group %in% group) %>%
    group_by(Group, Winged, Odor, PulseDuration) %>%
    summarise(N=n()) %>%
    group_by(Group, Winged) %>%
    summarise(N=unique(N))
  d_n <- unite(d_n, col = "Winged", sep = " (n = ")
  d_n$t <- ")"
  d_n <- unite(d_n, col = "Winged", sep = "")
  
  
  d_all <- d_mean
  
  
  data_se <- gather(d_se, time, se, t1:t5000, factor_key=TRUE)

  
  d_all$Odor <- factor(d_all$Odor)
  times <- times[times$Odor == "2-heptanone", ]
  
  data_long <- gather(d_all, time, Voltage, t1:t5000, factor_key=TRUE)
  data_long$time <- (as.numeric(gsub("t", "", data_long$time)) - 200)/1000
  data_long$SE <- data_se$se

  data_long$PulseDuration <- factor(data_long$PulseDuration, levels = c("50"), labels = c("10 Hz"))
  
  
  col <- c("grey", "#00DC00")
  
  data_segm1 <- data.frame(y=c(0, -0.02), yend=c(0.05, -0.02), x=c(Inf, 1), xend=c(Inf, 2),
                           PulseDuration = factor(150, levels = c(15, 30, 150)),
                           Group=factor(c("Six Mile", 'Mt Burns')),
                           Odor=c("2-heptanone"))
  label1 <- data.frame(x=c(Inf, 2), y=c(0.035, -0.035),
                       text=c("0.05 mV", '1 s'),
                       PulseDuration= factor(150, levels = c(15, 30, 150)),
                       Group=factor(c("Six Mile", 'Mt Burns')),
                       Odor=c("2-heptanone"))
  
  ## For Lug, we stimultaed with a different stimulus protocol were the stimulation
  ## with several pulses in a specific frequency was not 10 Hz and the total length
  ## o f stimulation was less than 3 seconds. We exclude these data here.
  dat <- data_long %>% filter(Group == creeks, Group != 'Lug')
  dat$Group <- factor(dat$Group, levels = c(creeks))
  
  W <- ggplot(dat, aes(x=time, y=Voltage)) + 
      geom_rect(data=times, inherit.aes=FALSE, fill = "grey95",
                aes(xmin=on, xmax=off, ymin=-Inf, ymax=Inf)) +
      # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
      geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
      geom_line(aes(colour = Winged)) +
      scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
      facet_grid(. ~ Group, scales = "free_y", drop = FALSE) +
      scale_y_continuous(breaks = pretty_breaks(n = 3)) +
      theme_bw() +
    labs(y = " ") +
      guides(fill = 'none') +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
      coord_cartesian(xlim =c(0, 3.5)) +
      theme(panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = "none",
            strip.background = element_rect(colour="NA", fill="NA"),
            strip.text.y = element_blank(),
            legend.title = element_blank(),
            axis.line = element_blank(),
            axis.ticks = element_blank(),
            axis.title.x = element_blank(),
            axis.title.y = element_text(size = 35),
            panel.border = element_blank(),
            axis.text=element_blank(),
            strip.text = element_blank(),
            text = element_text(size=18)) 
  W
  
  col <- c("blue")
  data_segm1 <- data.frame(y=0, yend=1, x=Inf, xend=Inf,
                           PulseDuration = factor(150, levels = c(15, 30, 150)),
                           Group=factor("Bee"),
                           Odor=c("2-heptanone"))
  label1 <- data.frame(x=Inf, y=0.5,
                       text=c("1 mV"),
                       PulseDuration= factor(150, levels = c(15, 30, 150)),
                       Group=factor("Bee"),
                       Odor=c("2-heptanone"))
  
  B <-  ggplot(data_long[data_long$Group %in% 'Bee', ], aes(x=time, y=Voltage)) + 
    geom_rect(data=times, inherit.aes=FALSE, fill = "grey95",
              aes(xmin=on, xmax=off, ymin=-Inf,ymax=Inf)) +
    # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
    geom_line(aes(colour = Winged)) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    facet_grid(Odor ~ Group, scales = "free_y") +
    scale_y_continuous(breaks = pretty_breaks(n = 3)) +
    theme_bw() +
    coord_cartesian(xlim =c(0, 3.5)) +
    # guides(fill = 'none') +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none",
          strip.background = element_rect(colour="NA", fill="NA"),
          panel.border = element_blank(),
          legend.title = element_blank(),
          axis.line = element_blank(),
          strip.text = element_blank(),
          axis.text=element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          text = element_blank()) 
  B
  
  
  data_segm1 <- data.frame(y=0, yend=1, x=3, xend=3,
                           PulseDuration = factor(150, levels = c(15, 30, 150)),
                           Group=factor("PID"),
                           Odor=c("2-heptanone"))
  label1 <- data.frame(x=3, y=0.5,
                       text=c("1 V"),
                       PulseDuration= factor(150, levels = c(15, 30, 150)),
                       Group=factor("PID"),
                       Odor=c("2-heptanone"))
  
  P <- ggplot(data_long[data_long$Group %in% 'PID', ], aes(x=time, y=Voltage)) + 
    geom_rect(data = times, inherit.aes = FALSE, fill = "grey95",
              aes(xmin = on, xmax=off, ymin = -Inf, ymax = Inf)) +
    # geom_ribbon(aes(ymin=Voltage + SE, ymax=Voltage - SE, fill=Winged), color=NA, alpha=0.3) +
    geom_hline(mapping = aes(yintercept = 0), color="black", size=0.3, linetype="dotted") +
    geom_line(aes(colour = Winged)) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    facet_grid(Odor ~ Group, scales = "free_y") +
    scale_y_continuous(breaks = pretty_breaks(n = 3)) +
    theme_bw() +
    coord_cartesian(xlim =c(0, 3.5)) +
    geom_segment(data=data_segm1,
                 aes(x=x, y=y, yend=yend, xend=xend)) +
    geom_text(data = label1, aes(x=x, y=y, label=text), size=geom.text.size, hjust=1.1) +
    # guides(fill = 'none') +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none",
          strip.background = element_rect(colour="NA", fill="NA"),
          panel.border = element_blank(),
          legend.title = element_blank(),
          axis.ticks = element_blank(),
          axis.text=element_blank(),
          axis.line = element_blank(),
          axis.title = element_blank(),
          strip.text = element_blank(),
          text = element_blank()) 
  P
  
  list(W = W, B = B, P = P)
  
}


FigS1C_new <- plot_freq(data2, c(creeks, 'Bee', 'PID'), Times)


## Power spectrum density analysis ####
## whole frequency band mean over groups of individuals ##
## Freq = 10 for individual ID ##
spec <- data.frame(x=as.numeric(NULL), y=as.numeric(NULL), sd=as.numeric(NULL), se=as.numeric(NULL),
                   Odor=as.character(NULL), Winged=as.character(NULL),
                   Group=as.character(NULL), TimeWindow=as.character(NULL))
spec_ind <- data.frame(y=as.numeric(NULL), 
                       yMean=as.numeric(NULL),
                       xMean=as.numeric(NULL),
                       zMean=as.numeric(NULL),
                       yMax=as.numeric(NULL),
                       xMax=as.numeric(NULL),
                       zMax=as.numeric(NULL),
                       id=as.numeric(NULL),
                       Odor=as.character(NULL), Winged=as.character(NULL),
                       Group=as.character(NULL), TimeWindow=as.character(NULL))


winged <- as.vector(t(unique(data2[, "Winged"])))

for(k in 1:length(winged)){
  group <- as.vector(t(unique(data2[data2$Winged==winged[k], "Group"])))
  
  for (l in 1:length(group)) {
    odors <- as.vector(t(unique(data2[data2$Winged==winged[k] & data2$Group==group[l], "Odor"])))
    
    for (i in 1:length(odors)){
      Hept <- data2 %>%
        filter(Odor==odors[i], Winged==winged[k], Group==group[l], (AntennaState=="alive" | is.na(AntennaState)))
      
      
      ### full length of 3 seconds ####
      spect <- apply(Hept[,paste0("t", 601:3600)],
                     MARGIN = 1,
                     FUN=function(x) spec.mtm(ts(x, frequency = 1000), taper = "sine", plot=F))
      
      spec <- rbind(spec, data.frame(x=spect[[1]]$freq[1:329],
                                     y=rowMeans(sapply(spect, '[[', 4))[1:329],
                                     sd=apply((sapply(spect, '[[', 4)),1, sd)[1:329],
                                     se=apply((sapply(spect, '[[', 4)),1, se)[1:329],
                                     Odor=odors[i],
                                     Winged=winged[k],
                                     Group=group[l],
                                     TimeWindow="T0-3"))
      
      # select frequencies between 9.5 and 10.5 and calculate the mean values
      sel_row_10 <- which(spect[[1]]$freq > 9 & spect[[1]]$freq < 11)
      sel_row_pre10 <- which(spect[[1]]$freq > 8 & spect[[1]]$freq < 9)
      sel_row_post10 <- which(spect[[1]]$freq > 11 & spect[[1]]$freq < 12)
      
      spec_ind <- rbind(spec_ind, data.frame(y=colSums(sapply(spect, '[[', 4)[sel_row_10,]),
                                             yMean=colMeans(sapply(spect, '[[', 4)[sel_row_10,]),
                                             xMean=colMeans(sapply(spect, '[[', 4)[sel_row_pre10,]),
                                             zMean=colMeans(sapply(spect, '[[', 4)[sel_row_post10,]),
                                             yMax=apply(sapply(spect, '[[', 4)[sel_row_10,], 2, max),
                                             xMax=apply(sapply(spect, '[[', 4)[sel_row_pre10,], 2, min),
                                             zMax=apply(sapply(spect, '[[', 4)[sel_row_post10,], 2, min),
                                             id=as.character(Hept$ID_new),
                                             Odor=odors[i],
                                             Winged=winged[k],
                                             Group=group[l],
                                             TimeWindow="T0-3"))
      
    }
  }
  
}



spec_sel <- spec %>%
  filter(Odor!="Blank", Odor!="Stonefly", Odor!="Propionic acid", TimeWindow=="T0-3")

spec_sel <- droplevels(spec_sel)

spec_sel$Odor <- factor(spec_sel$Odor, levels = c("2-heptanone",
                                                  "1-octanol",
                                                  "2-butanone"))


col <- c("#FF00FF", "#00DC00")



plot_freq <- function(data) {
  
  data <- data %>%
    filter(Odor=="2-heptanone")
  
  col <- c(fw_l, wr_l, "black")
  
  data$TimeWindow <- factor(data$TimeWindow, levels = levels(data$TimeWindow),
                            labels = "")
  
  odors <- data.frame("Odor"=levels(data$Odor))
  
  ## For Lug, we stimultaed with a different stimulus protocol were the stimulation
  ## with several pulses in a specific frequency was not 10 Hz and the total length
  ## o f stimulation was less than 3 seconds. We exclude these data here.
  dat <- data %>% filter(Group %in% creeks & Group != 'Lug')
  dat$Group <- factor(dat$Group, levels = c(creeks))
  
  W <- ggplot(dat,
              aes(x=x, y=log10(y))) + 
    geom_rect(mapping=aes(xmin=9.5, xmax=10.5,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    geom_ribbon(aes(ymin=log10(y + se), ymax=log10(y - se), fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    facet_grid(. ~ Group, drop = FALSE) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = col) +
    theme_bw() +
    scale_y_continuous(breaks = pretty_breaks(n=3)) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none",
          text = element_text(size=18),
          axis.title = element_blank(),
          axis.line = element_line(colour = "black"),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          axis.text=element_text(colour="black")) 
  W
  
  B <- ggplot(data %>% filter(Group %in% c("Bee")),
              aes(x=x, y=log10(y))) + 
    geom_rect(mapping=aes(xmin=9.5, xmax=10.5,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    geom_ribbon(aes(ymin=log10(y + se), ymax=log10(y - se), fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    facet_grid(. ~ Group) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = "blue") +
    theme_bw() +
    scale_y_continuous(breaks = pretty_breaks(n=3)) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none",
          text = element_text(size=18),
          axis.title = element_blank(),
          axis.line = element_line(colour = "black"),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          axis.text=element_text(colour="black"))
  B
  
  P <- ggplot(data %>% filter(Group %in% c("PID")),
              aes(x=x, y=log10(y))) + 
    geom_rect(mapping=aes(xmin=9.5, xmax=10.5,
                          ymin=-Inf, ymax=Inf),
              fill="grey90", alpha=1) +
    geom_ribbon(aes(ymin=log10(y + se), ymax=log10(y - se), fill=Winged),
                color=NA, alpha=0.3) +
    geom_line(aes(colour = Winged)) +
    facet_grid(. ~ Group) +
    scale_color_manual(aesthetics = c("colour", "fill"), values = "blue") +
    theme_bw() +
    scale_y_continuous(breaks = pretty_breaks(n=3)) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.position = "none",
          axis.title = element_blank(),
          text = element_text(size=18),
          axis.line = element_line(colour = "black"),
          strip.background = element_rect(colour="NA", fill="NA"),
          strip.text = element_blank(),
          panel.border = element_blank(),
          axis.text=element_text(colour="black"))
  P
  
  list1 <- list("traceW"=W,
                "traceB"=B,
                "traceP"=P)
  return(list1)
  
}

p5 <- plot_freq(data = spec_sel); p5
#create common x and y labels

FigS1D <- egg::ggarrange(p5$traceW, p5$traceB, p5$traceP,
                         bottom = textGrob("Frequency (Hz)", gp = gpar(fontsize=16), vjust = 0.5),
                         left = textGrob(expression(paste("Power (", mV^2/Hz, ") [log10]")),
                                         gp = gpar(fontsize=18), vjust = 0.5, rot = 90),
                         ncol = 3, nrow=1,
                         label.args = list(gp=gpar(fontface="bold", fontsize=18), hjust=0, vjust=1),
                         widths = c(5,1,1),
                         debug = FALSE)
FigS1D

FigS1C <- egg::ggarrange(FigS1C_new$W, FigS1C_new$B, FigS1C_new$P,
                         ncol = 3, nrow=1, heights = 1,
                         widths = c(5, 1, 1),
                         debug = FALSE)
FigS1C

Fig1AB <- cowplot::plot_grid(Fig1A, Fig1B, 
                             align = "hv", axis = "trbl",
                             labels = c("A", "B"),
                             rel_widths = c(1.2, 1.5),
                             nrow = 1, label_size = 18)

### Fig. 1 
cowplot::plot_grid(Fig1AB, Fig1C, nrow = 2, rel_heights = c(1, 1.3),
                   labels = c("", "C"), label_size = 18)


ggsave("Figure1.pdf", path = fig_path, width = 12, height = 10)


### Fig. 2
cowplot::plot_grid(Fig2A, Fig2B, nrow = 2,
                   labels = LETTERS[1:2], label_size = 18)


ggsave("Figure2.pdf", path = fig_path, width = 9, height = 12)


### Fig. S1
cowplot::plot_grid(FigS1A, FigS1B, NULL, FigS1C, FigS1D, nrow = 5, rel_heights = c(1.2, 1, 0.05, 0.3, 0.8),
                   labels = c('A', 'B', '', 'C', 'D'), label_size = 18)

ggsave("FigureS1.pdf", path = fig_path, width = 14, height = 14)