# load packages and data
library(tidyverse)
library(knitr)
# fitting psychometric function using quickpsy
library(quickpsy)
library(ez)
library(ggplot2)
library(ggsignif)
library(PairedData)
library(cowplot)
library(effsize)
library(dplyr)
library(gtools)
library(bayesplot)
library(gridExtra)
library(plyr)
library(ggpubr)
library(lmerTest)        # lmer function
library(car)             # contr.Sum function
library(sjPlot)          # tab_model function

theme_set(bayesplot::theme_default())

## define theme for whole report
mytheme <-  theme_bw()+theme(panel.border = element_blank(),
                              panel.grid.major = element_blank(), 
                              panel.grid.minor = element_blank(),
                              panel.background = element_blank(), 
                              axis.line = element_line(colour = "black"))+
  theme(strip.background = element_blank(),
        axis.text.x = element_text(color = "black", size = 10), 
        axis.text.y = element_text(color = "black", size = 10), 
        axis.title.x = element_text(color = "black", size = 11, vjust = -1),
        axis.title.y = element_text(color = "black", size = 11),
        plot.title = element_text(color = "black", size = 11),
        legend.text=element_text(size=9),axis.ticks.length=unit(.10, "cm"))

colors_plot <- c("#185C9B","#D81E3D")

##########################################################################
# load behavior data
dat_EEG_exp2 = read.csv('data/allData_exp2_beh.csv') 
dat_EEG_exp1 = read.csv('data/allData_exp1_beh.csv')
eeg_dat_all <- data.frame(rbind(dat_EEG_exp1, dat_EEG_exp2))

# load ERP data
allAverageDat_exp1 <- read.csv(paste0(getwd(), '/data/allAverageDat_cnv_exp1.csv'))
allAverageDat_exp2 <- read.csv(paste0(getwd(), '/data/allAverageDat_cnv_exp2.csv'))

data_ica_p2_exp1 <- read.csv(paste0(getwd(), '/data/allAverageDat_pos_exp1.csv'))
data_ica_pos_exp2 <- read.csv(paste0(getwd(), '/data/allAverageDat_pos_exp2.csv'))

# load peak data
exp1_negPeak_lat <- read.csv(paste0(getwd(), '/data/peak_negative_cnv_exp1.csv'))
exp2_negPeak_lat <- read.csv(paste0(getwd(), '/data/peak_negative_cnv_exp2.csv'))

exp1_pos_peaks <- read.csv(paste0(getwd(), '/data/peak_positive_offset_exp1.csv'))
exp2_pos_peaks <- read.csv(paste0(getwd(), '/data/peak_positive_offset_exp2.csv'))

# CNV end-turning point data
turn_exp1= read.csv("data/width_CNV_exp1.csv")
turn_exp2= read.csv("data/width_CNV_exp2.csv")

## calulate geometric mean
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}