# Need to visually inspect more the data, see what to do with outliers both at participant and trial level. This is because it influence the inverse efficiency analysis
####### import packages ######
# load packages and functions
set.seed(1000)
library(lme4)
library(lmerTest)
library(dplyr)
library(ggplot2)
sem_fun = function(x) sd(x)/sqrt(length(x)) # standard error of the mean
library(tidyr)
library(dplyr)
library("reshape2")
library(rstatix)
library(readr)
library(viridis)
library(ggridges)
plot_pp_outliers = function(...){                                     
  # Participants median and accuracy plotted with interquartile range
  ciao = dataset %>% 
    group_by(subj_nr) %>%
    summarise(acc = mean(correct))
  
  ciao2 = dataset[as.logical(dataset$correct),] %>% 
    group_by(subj_nr) %>%
    summarise(rt = median(rt))
  
  ciao$rt = ciao2$rt
  ciao = ciao %>%
    ungroup()%>%
    mutate(upper_bound_rt = quantile(rt,0.75) + IQR(rt)*1.5, rt,  lower_bound_acc = quantile(acc,0.25) - IQR(acc)*1.5, acc, outlier = acc < lower_bound_acc | rt > upper_bound_rt,n_pp = sum(!outlier))
  
  
  
  p1 <- ggplot(ciao, aes(acc, rt, col = as.factor(outlier))) + geom_point() + theme_bw() + scale_x_reverse() +
    geom_hline(yintercept = ciao$upper_bound_rt[1],linetype='dashed', col = 'red',size=1) +
    geom_vline(xintercept = ciao$lower_bound_acc[1],linetype='dashed', col = 'red',size=1)+
    labs(title = 'Reation Times and Accuracy per participant',
         subtitle = paste0('Lines indicate outlier bounds\nNumber of participants remaining: ',ciao$n_pp))+
    theme_grey()+
    theme(legend.position = "none")+
    geom_text(aes(label=as.character(subj_nr)),hjust=0,vjust=0)+
    xlab("Accuracy (flipped)") + ylab("Median reaction times")
  
  
  return(ggExtra::ggMarginal(
    p = p1,
    type = 'boxplot',
    margins = 'both',
    size = 5,
    colour = 'gray4',
    fill = '#F50A0AB8'
  ) )
} 
source("IQR_custom.R") # function for Inter-Quartile Range
library(dplyr)
library(ggplot2)
library(effsize)
sem_fun = function(x) sd(x)/sqrt(length(x)) # standard error of the mean
require(dplyrExtras)
# modify violin
{
  "%||%" <- function(a, b) {
    if (!is.null(a)) a else b
  }
  
  geom_flat_violin <- function(mapping = NULL, data = NULL, stat = "ydensity",
                               position = "dodge", trim = TRUE, scale = "area",
                               show.legend = NA, inherit.aes = TRUE, ...) {
    layer(
      data = data,
      mapping = mapping,
      stat = stat,
      geom = GeomFlatViolin,
      position = position,
      show.legend = show.legend,
      inherit.aes = inherit.aes,
      params = list(
        trim = trim,
        scale = scale,
        ...
      )
    )
  }
  
  GeomFlatViolin <-
    ggproto("GeomFlatViolin", Geom,
            setup_data = function(data, params) {
              data$width <- data$width %||%
                params$width %||% (resolution(data$x, FALSE) * 0.9)
              
              # ymin, ymax, xmin, and xmax define the bounding rectangle for each group
              data %>%
                group_by(group) %>%
                mutate(ymin = min(y),
                       ymax = max(y),
                       xmin = x - width,
                       xmax = x)
            },
            
            draw_group = function(data, panel_scales, coord) {
              
              # # define what should be left and what right 
              data$side_plot[data$group %% 2 == 0] = -1 #right
              data$side_plot[data$group %% 2 != 0] = 1 #left
              
              # Find the points for the line to go all the way around
              data <- data %>%
                group_by(group) %>%
                mutate(xmaxv = x,
                       xminv = x + (side_plot * violinwidth) * (xmin - x))
              
              
              
              
              
              # Make sure it's sorted properly to draw the outline
              newdata <- rbind(plyr::arrange(transform(data, x = xminv), y),
                               plyr::arrange(transform(data, x = xmaxv), -y))
              
              # Close the polygon: set first and last point the same
              # Needed for coord_polar and such
              newdata <- rbind(newdata, newdata[1,])
              
              ggplot2:::ggname("geom_flat_violin", GeomPolygon$draw_panel(newdata, panel_scales, coord))
            },
            
            draw_key = draw_key_polygon,
            
            default_aes = aes(weight = 1, colour = "grey20", fill = "white", size = 0.5,
                              alpha = NA, linetype = "solid"),
            
            required_aes = c("x", "y")
    )
}

######## Load_data ########
dataset <- read_csv("exp2_dataset.csv")

number_of_subjs = length(unique(dataset$PROLIFIC_PID))

# check number of trials
if (length(dataset$PROLIFIC_PID)/length(unique(dataset$PROLIFIC_PID)) == 540){
  
  n_trials_per_subj = 540
  
} else { print("ERROR: number of trials are not 540 per subject")}

# assign number to participants 
dataset$subj_nr =  rep(1:number_of_subjs, each = n_trials_per_subj)
dataset$subj_nr <- as.factor(dataset$subj_nr)

####### LAG ######### 


# actual lag inclusion
{
  # check for keyboard response lagging (time from mask to accepted response, shoud be 0)
  summary(dataset$time_keyboard_response - dataset$time_mask) # mask duration should be 0
  
  # check for duration of target and delay (expected 60 + 100 ms)
  summary(dataset$time_delay_target - dataset$time_target) # duration of target presentation should be 60
  summary(dataset$time_mask - dataset$time_delay_target) # duration of delay should be 100
  
  
  # DECISION #
  # I will add to the reaction times both the time of target presentation and the time of the delay
  
  dataset$target_duration = dataset$time_delay_target - dataset$time_target
  dataset$delay_duration =  dataset$time_mask - dataset$time_delay_target
  dataset$mask_duration = dataset$time_keyboard_response - dataset$time_mask
  
  # Actual reaction times
  dataset$reaction_time =  dataset$response_time_keyboard_response + dataset$delay_duration + dataset$target_duration + dataset$mask_duration
  
  summary(dataset$reaction_time)
  
  dataset$rt = dataset$reaction_time
}


# cut <1700
dataset = dataset[dataset$rt < 1700,]

# Identify outliers
windows()
plot_pp_outliers()

# Exclude outliers
dataset = dataset[!dataset$subj_nr %in% c(13,16,32,19,10),]

dataset$subj_nr = factor(dataset$subj_nr)
dataset$block_number = as.factor(dataset$block_number)

windows()
plot_pp_outliers()

######## prepare conditions #########

valid_trials = dataset$target_pos == "valid"
sum(valid_trials)/number_of_subjs

# tranform coordinates of cue pos in side (left, right)
dataset$cue_pos_x
dataset$cue_side[dataset$cue_pos_x < 0] = 'left'
dataset$cue_side[dataset$cue_pos_x > 0] = 'right'


######## Object based effect and reward ###########
obj_eff_rwd_data = dataset[dataset$target_pos != "valid",] # taking into account only invalid trials

obj_eff_rwd_data$target_pos = as.factor(as.character(obj_eff_rwd_data$target_pos)) #### added to trim valid level

obj_eff_rwd_data$target_side = obj_eff_rwd_data$target_x == 192
obj_eff_rwd_data$target_side[obj_eff_rwd_data$target_side] = "right"
obj_eff_rwd_data$target_side[obj_eff_rwd_data$target_side == "FALSE"] = "left"

obj_eff_rwd_data$rew_trg_rel = rep("placeholder",length(obj_eff_rwd_data$PROLIFIC_PID))

obj_eff_rwd_data$rew_trg_rel[obj_eff_rwd_data$side_high_reward == "to_be_assigned"] = "low_reward"

obj_eff_rwd_data$rew_trg_rel[(obj_eff_rwd_data$target_side == obj_eff_rwd_data$side_high_reward)] = "congruent"
obj_eff_rwd_data$rew_trg_rel[(obj_eff_rwd_data$target_side != obj_eff_rwd_data$side_high_reward & obj_eff_rwd_data$side_high_reward != "to_be_assigned")] = "incongruent"

obj_eff_rwd_data$rew_trg_rel = as.factor(obj_eff_rwd_data$rew_trg_rel)
obj_eff_rwd_data$target_pos = as.factor(obj_eff_rwd_data$target_pos)

obj_eff_rwd_data = obj_eff_rwd_data %>%
  mutate(rew_trg_rel = forcats::fct_relevel(rew_trg_rel, 
                                            "low_reward", "congruent", "incongruent"))


obj_eff_rwd_data = obj_eff_rwd_data %>% rename(subj=subj_nr, rwd = rew_trg_rel, trg = target_pos)

# Accuracy per participant
acc = dataset %>%
  group_by(PROLIFIC_PID)%>%
  summarise(accuracy = mean(correct))
summary(acc$accuracy)

# Number of trials in the smallest cell per subj
obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
  group_by(rwd, trg,subj) %>%
  summarise(rt = mean(rt), n=n()) %>%
  group_by(rwd, trg) %>%
  summarise(rt = mean(rt), n=mean(n)) %>% print(n = 100)


######## plots ##########

## new plot ####

{
  plot_data = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt), n=n())
  
  plot_data = transform(plot_data, drwd = ifelse(trg == "diff_obj", 
                                                 as.numeric(rwd)-.125, # 0.25 with position_dodge = 1
                                                 as.numeric(rwd)+.125 ) )
  plot_data2 = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(rwd, trg) %>%
    summarise(rt = mean(rt), n=n())
  
  plot_data2 = transform(plot_data2, drwd = ifelse(trg == "diff_obj", 
                                                   as.numeric(rwd)-.125,
                                                   as.numeric(rwd)+.125 ) )
  
  
  windows()
  p = ggplot(data=plot_data, aes(x=rwd, y=rt, fill=trg))+
    geom_flat_violin(position = position_dodge(width=1.5))+
    geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), position = position_dodge(width=1) )+
    geom_point(aes(group=interaction(trg, rwd), shape=trg),alpha= 3/10, position = position_dodge(width=0.5))+
    geom_line(aes(x=drwd, y=rt, group=interaction(subj, rwd)),alpha= 3/10)+
    geom_point(data=plot_data2, aes(group=interaction(trg, rwd), shape=trg), size=3, position = position_dodge(width=0.5))+
    geom_line(data=plot_data2, aes(x=drwd, y=rt, group=rwd),size = 1)+
    scale_shape_manual(values=c(15,16))+
    theme_minimal()+
    theme(plot.title = element_text(hjust = 0.5))+
    coord_cartesian(ylim = c(500,1200))+
    scale_y_continuous(breaks=seq(500, 1200, 100))
  p
} 



# plot RT, rew-trg and cue-trg ####
{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, trg) %>%
    group_by(rwd, trg) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, trg)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=rt, fill=trg))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    geom_errorbar(aes(ymin=rt-isv, ymax=rt+isv), width=.2,position=position_dodge(.5))+
    coord_cartesian(ylim = c(600,950))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    scale_y_continuous(breaks=seq(600, 960, 20))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}

obj_eff_rwd_data$block_number = factor(obj_eff_rwd_data$block_number)
# block number - line plot ####
{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj,block_number) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj, block_number,rwd) %>% 
    summarize(rt = rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd,block_number) %>%
    group_by(rwd, block_number) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, block_number)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot( aes(x=as.numeric(block_number), y=rt, colour=rwd)) + geom_point() + geom_line()+
    geom_ribbon(aes(ymin=rt-isv, ymax=rt+isv), linetype=2, alpha=0.1)+
    theme_minimal()+
    #    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    #    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}

# block number and subj number - line plot ####
{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj,block_number) %>%
    summarise(rt = mean(rt)) %>%
    group_by(rwd, block_number, subj)%>%
    summarize(rt =  rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>%
    ggplot( aes(x=as.numeric(block_number), y=rt, colour=subj)) + geom_point() + geom_line()+
    theme_minimal()+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")+
    facet_grid(. ~ rwd)
  p
}


###### inverse efficiency score ###########
{
  
  # inverse efficiency function
  inv_efficiency = function(...){
    obj_eff_rwd_data = obj_eff_rwd_data %>% # first calculate accuracy per condition per subj
      group_by(rwd, subj, trg)%>%
      mutate(accuracy = mean(correct))%>% # lenght remain the same
      as.data.frame()
    
    
    obj_eff_rwd_data_correct = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] # keep only correct trials
    
    obj_eff_rwd_data_meanSj = obj_eff_rwd_data_correct %>% # mean by rwd, subj 
      group_by(rwd, subj, trg)%>%
      summarise(rt = mean(rt), accuracy=mean(accuracy))%>% #mean accuracy doesn't affect the value because it's already a measure by reward cue rel, subj and block type
      as.data.frame()
    
    obj_eff_rwd_data_meanSj = obj_eff_rwd_data_meanSj %>% #  create inverse efficiency score
      group_by(rwd, subj, trg)%>%
      mutate(inv_eff = rt/accuracy)%>%
      as.data.frame()
    
    return(obj_eff_rwd_data_meanSj)
  }
  obj_eff_rwd_data_meanSj = inv_efficiency() 
  
  
}

# I do not compute outliers on inv efficiency 
# plot inverse efficiency score, rew-trg and cue-trg ######
{
  windows()
  p = obj_eff_rwd_data_meanSj %>%
    group_by(rwd, trg)%>%
    summarise(stand_dev = sd(inv_eff), inv_eff = mean(inv_eff))  %>% #  stand_dev = sd(inv_eff), put it before inv_eff otherwise won't work
    ggplot(aes(x=rwd, y=inv_eff, fill=trg))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    #    geom_errorbar(aes(ymin=inv_eff-stand_dev, ymax=inv_eff+stand_dev), width=.2,position=position_dodge(.5))+
    #    coord_cartesian(ylim = c(500,1800))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    ggtitle("Object-based effect\n(error bar: standard deviation)")+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("Inverse Efficiency Score")
  p
}




# plot orientation ####
{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, trg) %>%
    group_by(rwd, trg) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, trg)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=rt, fill=trg))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    geom_errorbar(aes(ymin=rt-isv, ymax=rt+isv), width=.2,position=position_dodge(.5))+
    coord_cartesian(ylim = c(600,950))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    scale_y_continuous(breaks=seq(600, 960, 20))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}


length(unique(obj_eff_rwd_data_meanSj$subj))
length(unique(dataset$PROLIFIC_PID))


# plot accuracy, rew-trg and cue-trg #####
{
  windows()
  p = obj_eff_rwd_data_meanSj %>%
    group_by(rwd, trg)%>%
    summarise(accuracy = mean(accuracy))  %>% #  stand_dev = sd(accuracy), put it before accuracy otherwise won't work
    ggplot(aes(x=rwd, y=accuracy, fill=trg))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    # geom_errorbar(aes(ymin=accuracy-stand_dev, ymax=accuracy+stand_dev), width=.2,position=position_dodge(.5))+
    # coord_cartesian(ylim = c(500,700))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    ggtitle("Object-based effect\n(error bar: intra-subj variability)")+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("Accuracy")
  p
}


############################ STATS ###################################
#### 2 paired samples t.test: cue-trg --- low reward ####

### only LOW REW
obj_eff_low_rew = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward == "to_be_assigned",]
obj_eff_low_rew$rwd = factor(obj_eff_low_rew$rwd)  # relevel


per_anova_earn = obj_eff_low_rew[as.logical(obj_eff_low_rew$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)

# RT
A_ttest_same = per_anova_earn$rt[per_anova_earn$trg == "same_obj"]
A_ttest_diff = per_anova_earn$rt[per_anova_earn$trg == "diff_obj"]

t.test(A_ttest_same, A_ttest_diff, paired = T) # t test
ttestBF(A_ttest_same, A_ttest_diff, paired = T) # BF

mean(A_ttest_same)
mean(A_ttest_diff)

library(effsize)
cohen.d(A_ttest_same, A_ttest_diff, paired = T)

# inv efficiency
inv_same = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "same_obj" & obj_eff_rwd_data_meanSj$rwd == "low_reward"]
inv_diff = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "diff_obj" & obj_eff_rwd_data_meanSj$rwd == "low_reward"]

t.test(inv_same, inv_diff, paired = T) # t test
ttestBF(inv_same, inv_diff, paired = T) # BF

mean(inv_same)
mean(inv_diff)

## ACC ##

supporto = obj_eff_low_rew %>%
  group_by(subj, rwd, trg) %>%
  summarise(acc = mean(correct))

obj_eff_low_rew %>%
  group_by(subj, rwd, trg) %>%
  summarise(acc = mean(correct)) %>%
  group_by(rwd, trg) %>%
  summarise(acc = mean(acc))

supporto = ungroup(supporto)
str(supporto)

library(MASS)

b <- boxcox(lm(supporto$acc ~ 1)) # boxcox on ACC
lambda <- b$x[which.max(b$y)] # Exact lambda
lambda
supporto$correct_ed <- (supporto$acc ^ lambda - 1) / lambda
hist(supporto$correct_ed)
hist(supporto$acc)

t.test(supporto$correct_ed[supporto$trg == 'same_obj'],supporto$correct_ed[supporto$trg == 'diff_obj'], paired = T)

cohen.d(supporto$correct_ed[supporto$trg == 'same_obj'],supporto$correct_ed[supporto$trg == 'diff_obj'], paired = T)



### ONLY congruent
obj_eff_cong_rwd = obj_eff_rwd_data[obj_eff_rwd_data$rwd == "congruent",]
obj_eff_cong_rwd$rwd = factor(obj_eff_cong_rwd$rwd)  # relevel 


per_anova_earn = obj_eff_cong_rwd[as.logical(obj_eff_cong_rwd$correct) ,] %>%
  group_by(rwd, subj, trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)


A_ttest_same = per_anova_earn$rt[per_anova_earn$trg == "same_obj"]
A_ttest_diff = per_anova_earn$rt[per_anova_earn$trg == "diff_obj"]

t.test(A_ttest_same, A_ttest_diff, paired = T) # t test
ttestBF(A_ttest_same, A_ttest_diff, paired = T) # BF


# inv efficiency
cong_same = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "same_obj" & obj_eff_rwd_data_meanSj$rwd == "congruent"]
cong_diff = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "diff_obj" & obj_eff_rwd_data_meanSj$rwd == "congruent"]

t.test(cong_same, cong_diff, paired = T) # t test
b = ttestBF(cong_same, cong_diff, paired = T) # BF
1/b

### ONLY incongruent

obj_eff_incong_rwd = obj_eff_rwd_data[obj_eff_rwd_data$rwd == "incongruent",]
obj_eff_incong_rwd$rwd = factor(obj_eff_incong_rwd$rwd)  # relevel 

per_anova_earn = obj_eff_incong_rwd[as.logical(obj_eff_incong_rwd$correct),] %>%
  group_by(rwd, subj, trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)


A_ttest_same = per_anova_earn$rt[per_anova_earn$trg == "same_obj"]
A_ttest_diff = per_anova_earn$rt[per_anova_earn$trg == "diff_obj"]

t.test(A_ttest_same, A_ttest_diff, paired = T) # t test
ttestBF(A_ttest_same, A_ttest_diff, paired = T) # BF

# inv efficiency
incong_same = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "same_obj" & obj_eff_rwd_data_meanSj$rwd == "incongruent"]
incong_diff = obj_eff_rwd_data_meanSj$inv_eff[obj_eff_rwd_data_meanSj$trg == "diff_obj" & obj_eff_rwd_data_meanSj$rwd == "incongruent"]

t.test(incong_same, incong_diff, paired = T)
ttestBF(A_ttest_same, A_ttest_diff, paired = T) # BF

#### 2X2 anova: congruent(rew-trg), cue-trg --- high rew vs low rew baseline ####

# congruent vs baseline

obj_eff_cong_rwd = obj_eff_rwd_data[obj_eff_rwd_data$rwd != "incongruent",]
obj_eff_cong_rwd$rwd = factor(obj_eff_cong_rwd$rwd)  # relevel 


per_anova_earn = obj_eff_cong_rwd[as.logical(obj_eff_cong_rwd$correct) ,] %>%
  group_by(subj, rwd,trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)

res.aov <- anova_test(data = per_anova_earn, dv = rt, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 10000)
windows()
plot(bf)
title("2X2\ncongruent vs baseline\nanovaBF", adj = 0)
against_best = bf/bf
plot(against_best[,3])
title("2X2\ncongruent vs baseline\nanovaBF", adj = 0)
bf


# inv efficiency

per_anova_earn = obj_eff_rwd_data_meanSj[obj_eff_rwd_data_meanSj$rwd != "incongruent",]
per_anova_earn$rwd = factor(per_anova_earn$rwd)  # relevel 


res.aov <- anova_test(data = per_anova_earn, dv = inv_eff, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

#### 2X2 anova: incongruent(rew-trg), cue-trg --- high rew vs low rew baseline ####

# incongruent vs baseline

obj_eff_incong_rwd = obj_eff_rwd_data[obj_eff_rwd_data$rwd != "congruent",]
obj_eff_incong_rwd$rwd = factor(obj_eff_incong_rwd$rwd)  # relevel


per_anova_earn = obj_eff_incong_rwd[as.logical(obj_eff_incong_rwd$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)

res.aov <- anova_test(data = per_anova_earn, dv = rt, wid = c(subj),  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)
bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 1000000)
windows()
plot(bf)
title("2X2\nincongruent vs baseline\nanovaBF", adj = 0)

against_best = bf/bf
windows()
plot(against_best[4,3])
title("2X2\nincongruent vs baseline\nanovaBF", adj = 0)

# model_1 = lmer(rt ~ rwd * trg + (1 + block_number|subj),
#   data = obj_eff_incong_rwd)
# 
# summary(model_1)

# inv efficiency
per_anova_earn = obj_eff_rwd_data_meanSj[obj_eff_rwd_data_meanSj$rwd != "congruent",]
per_anova_earn$rwd = factor(per_anova_earn$rwd) # relevel 


res.aov <- anova_test(data = per_anova_earn, dv = inv_eff, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

#### 2X2 anova: rew-trg and cue-trg --- high reward #####

# without low reward baseline

obj_eff_high_rwd = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward != "to_be_assigned",]
obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd)  # relevel with only high rwd trials

## Statistical Model

per_anova_earn = obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)

res.aov <- anova_test(data = per_anova_earn, dv = rt, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 1000000)
windows()
plot(bf)
title("2X2\nHigh reward\nanovaBF", adj = 0)

against_best = bf/bf
windows()
plot(against_best[4,3])
title("2X2\nincongruent vs baseline\nanovaBF", adj = 0)


# inv efficiency
per_anova_earn = obj_eff_rwd_data_meanSj[obj_eff_rwd_data_meanSj$rwd != "low_reward",]
per_anova_earn$rwd = factor(per_anova_earn$rwd)  # relevel 


res.aov <- anova_test(data = per_anova_earn, dv = inv_eff, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)


##### ACC #####


obj_eff_high_rwd %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct)) %>%
  group_by(rwd,trg)%>%
  summarise(acc = mean(acc))

per_anova_earn = obj_eff_high_rwd %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct))

per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)


library(MASS)

b <- boxcox(lm(per_anova_earn$acc ~ 1)) # boxcox on ACC
lambda <- b$x[which.max(b$y)] # Exact lambda
lambda
per_anova_earn$correct_ed <- (per_anova_earn$acc ^ lambda - 1) / lambda
hist(per_anova_earn$correct_ed)
hist(per_anova_earn$acc)




res.aov <- anova_test(data = per_anova_earn, dv = correct_ed, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

#main rwd
obj_eff_high_rwd %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct)) %>%
  group_by(rwd)%>%
  summarise(acc = mean(acc))

#main trg
obj_eff_high_rwd %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct)) %>%
  group_by(trg)%>%
  summarise(acc = mean(acc))

# interaction
obj_eff_high_rwd %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct)) %>%
  group_by(trg, rwd)%>%
  summarise(acc = mean(acc))



#### 3X2 ####
obj_eff_rwd_data

per_anova_earn = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(rt = mean(rt))

#str(per_anova_earn)
per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)

res.aov <- anova_test(data = per_anova_earn, dv = rt, wid = subj,  within = c(rwd, trg), effect.size = "pes")
get_anova_table(res.aov)

bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 1000000)
bf
windows()
plot(bf)
title("3X2\nanovaBF", adj = 0)

against_best = bf/bf
against_best

windows()
plot(against_best[4,3])
title("2X2\nincongruent vs baseline\nanovaBF", adj = 0)