# libs and funs
# 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)
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'
  ) )
}                                           
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 Jatos results
dataset = read_csv("D:/Motivated Cognition Lab/UBIRA_JEP_paper/exp1/exp1_dataset.csv")

# Delete practice
dataset = dataset[dataset$block_number != 0,]

# assign number to participants (trials has to be 648)
number_of_subjects = length(unique(dataset$PROLIFIC_PID))
dataset$number =  rep(1:number_of_subjects, each = 648)
dataset$number <- as.factor(dataset$number)
dataset$subj_nr= dataset$number

levels(dataset$number)
dataset$count_trial_sequence = rep(1:648,length(levels(dataset$subj_nr)))

{
  # check for keyboard response lagging (time from mask to accepted response, shoud be 0)
  # to do so we merge the between condition variables in a single variable
  dataset$time_keyboard_response[is.na(dataset$time_keyboard_response)] = dataset$time_keyboard_response_horiz[!is.na(dataset$time_keyboard_response_horiz)]
  dataset$time_mask[is.na(dataset$time_mask)] = dataset$time_mask_horiz[!is.na(dataset$time_mask_horiz)]
  
  
  summary(dataset$time_keyboard_response - dataset$time_mask) # mask duration should be 0
  
  # check for duration of target and delay (expected 60 + 100 ms)
  dataset$time_delay_target[is.na(dataset$time_delay_target)] = dataset$time_delay_target_horiz[!is.na(dataset$time_delay_target_horiz)]
  dataset$time_target[is.na(dataset$time_target)] = dataset$time_target_horiz[!is.na(dataset$time_target_horiz)]
  
  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
  
  # I will add to the reaction times both the time of target presentation and the time of the delay
  
  dataset$mask_duration = dataset$time_keyboard_response - dataset$time_mask
  dataset$target_duration = dataset$time_delay_target - dataset$time_target
  dataset$delay_duration =  dataset$time_mask - dataset$time_delay_target
  
  # is there any NAs?
  any(is.na(dataset$mask_duration))
  any(is.na(dataset$target_duration))
  any(is.na(dataset$delay_duration))
  
  # where is the NA
  unique(dataset$PROLIFIC_PID[is.na(dataset$target_duration)])
  
  # merge the between condition response time in a single variable
  dataset$response_time_keyboard_response[is.na(dataset$response_time_keyboard_response)] = dataset$response_time_keyboard_response_horiz[!is.na(dataset$response_time_keyboard_response_horiz)]
  
  # 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
  
}

# nothing over 1700 ms
dataset$rt_outliers_logic = dataset$rt<1700

dataset %>% 
  group_by(subj_nr) %>%
  summarise(outl = (1-sum(rt_outliers_logic)/length(rt_outliers_logic))*100) %>%
  group_by()%>%
  summarise(sd = sd(outl), outl = mean(outl))

# outier => deleted
dataset = dataset[dataset$rt_outliers_logic,]

windows()
plot_pp_outliers()


# 60, 61, 59 are outliers (use plot_pp_outliers())
dataset = dataset[!dataset$subj_nr %in% c(60,61,59),]
number_of_subjects = length(unique(dataset$PROLIFIC_PID))
dataset$PROLIFIC_PID = factor(dataset$PROLIFIC_PID)
dataset$number = factor(dataset$number)

# define Valid VS Invalid
dataset$space = ifelse(dataset$target_pos == "valid", 'valid', 'invalid')  
dataset$rectangle_orientation = factor(dataset$rectangle_orientation)
dataset$block_number = factor(dataset$block_number)


###### RT ######
per_ttest = dataset[dataset$accuracy == 100,] %>%
  subset(target_pos != "valid") %>%
  mutate(factor(target_pos)) %>%
  group_by(target_pos,subj_nr)%>%
  summarise(rt = mean(reaction_time))

t.test(per_ttest$rt[per_ttest$target_pos == "same_obj"], per_ttest$rt[per_ttest$target_pos == "diff_obj"], paired = T)
mean(per_ttest$rt[per_ttest$target_pos == "same_obj"]-per_ttest$rt[per_ttest$target_pos == "diff_obj"])
sd(per_ttest$rt[per_ttest$target_pos == "same_obj"]-per_ttest$rt[per_ttest$target_pos == "diff_obj"])
cohen.d(per_ttest$rt[per_ttest$target_pos == "same_obj"], per_ttest$rt[per_ttest$target_pos == "diff_obj"], paired = T)

mean(per_ttest$rt[per_ttest$target_pos == "same_obj"])
mean(per_ttest$rt[per_ttest$target_pos == "diff_obj"])


##### violin plot ######


{
  plot_data = dataset[as.logical(dataset$correct),] %>%
    group_by(subj_nr, target_pos) %>%
    summarise(rt = mean(rt), n=n())
  
  
  plot_data2 = dataset[as.logical(dataset$correct),] %>%
    group_by(subj_nr, target_pos) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(target_pos) %>%
    summarise(rt = mean(rt), n=n())
  
  windows()
  p = ggplot(data=plot_data, aes(x=target_pos, y=rt, fill=target_pos))+
    geom_violin(draw_quantiles = c(0.25, 0.5, 0.75))+
    geom_point(aes(group=target_pos),alpha= 3/10)+
    geom_line(aes(x=target_pos, y=rt, group=subj_nr),alpha= 3/10)+
    geom_point(data=plot_data2, aes(group=target_pos), size=3)+
    geom_line(data=plot_data2, aes(x=target_pos, y=rt, group=target_pos),size = 1)+
    scale_shape_manual(values=c(15))+
    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 space- and object-based effect 
{
  windows()
  p = dataset[as.logical(dataset$correct),] %>%
    group_by(subj_nr, target_pos) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj_nr) %>%
    mutate(subj_mean = mean(rt), target_pos) %>%
    group_by(target_pos) %>%
    mutate(group_mean = mean(rt), subj_nr) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(target_pos)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=target_pos, y=rt, fill=target_pos))+
    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(625,700))+
    theme_minimal()+
    scale_x_discrete(labels= c("different", "same", "valid"))+
    scale_fill_discrete(name = "Cue & Traget", labels = c("different", "same", "valid"))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(factor(dataset$subj_nr)))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Cue & Traget") + ylab("RT")
  p
}

# space
supporto = dataset[as.logical(dataset$correct),] %>%
  group_by(subj_nr, space) %>%
  summarise(rt = mean(rt))

dataset[as.logical(dataset$correct),] %>%
  group_by(subj_nr, space) %>%
  summarise(rt = mean(rt))%>%
  group_by(space) %>%
  summarise(rt = mean(rt))

t.test(supporto$rt[supporto$space == 'valid'],supporto$rt[supporto$space == 'invalid'], paired = T)
cohen.d(supporto$rt[supporto$space == 'valid'],supporto$rt[supporto$space == 'invalid'], paired = T)

# object
supporto = dataset[as.logical(dataset$correct),] %>%
  subset(target_pos != 'valid') %>%
  group_by(subj_nr, target_pos) %>%
  summarise(rt = mean(rt))
t.test(supporto$rt[supporto$target_pos == 'same_obj'],supporto$rt[supporto$target_pos == 'diff_obj'], paired = T)

dataset[as.logical(dataset$correct),] %>%
  subset(target_pos != 'valid') %>%
  group_by(subj_nr, target_pos) %>%
  summarise(rt = mean(rt)) %>%
  group_by(target_pos) %>%
  summarise(rt = mean(rt))

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


###### ACC ######
# plot Accuracy space- and object-based effect 

{
  windows()
  p = dataset %>%
    group_by(subj_nr, target_pos) %>%
    summarise(acc = mean(acc)) %>%
    group_by(subj_nr) %>%
    mutate(subj_mean = mean(acc), target_pos) %>%
    group_by(target_pos) %>%
    mutate(group_mean = mean(acc), subj_nr) %>%
    group_by() %>%
    mutate(isv = (acc - subj_mean + group_mean)) %>%
    group_by(target_pos)%>%
    summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
    ggplot(aes(x=target_pos, y=acc, fill=target_pos))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    geom_errorbar(aes(ymin=acc-isv, ymax=acc+isv), width=.2,position=position_dodge(.5))+
    coord_cartesian(ylim = c(70,100))+
    theme_minimal()+
    scale_x_discrete(labels= c("different", "same", "valid"))+
    scale_fill_discrete(name = "Cue & Traget", labels = c("different", "same", "valid"))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(factor(dataset$subj_nr)))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Cue & Traget") + ylab("RT")
  p
}

# validity
supporto = dataset %>%
  group_by(subj_nr, space) %>%
  summarise(acc = mean(acc))
# t.test(supporto$acc[supporto$space == 'valid'],supporto$acc[supporto$space == 'invalid'], paired = T)
# cohen.d(supporto$acc[supporto$space == 'valid'],supporto$acc[supporto$space == 'invalid'], paired = T)

dataset %>%
  group_by(subj_nr, space) %>%
  summarise(acc = mean(acc)) %>%
  group_by(space) %>%
  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$space == 'valid'],supporto$correct_ed[supporto$space == 'invalid'], paired = T)
cohen.d(supporto$correct_ed[supporto$space == 'valid'],supporto$correct_ed[supporto$space == 'invalid'], paired = T)


# object based effect

supporto = dataset %>%
  subset(target_pos != 'valid') %>%
  group_by(subj_nr, target_pos) %>%
  summarise(acc = mean(acc))

dataset %>%
  subset(target_pos != 'valid') %>%
  group_by(subj_nr, target_pos) %>%
  summarise(acc = mean(acc)) %>%
  group_by(target_pos) %>%
  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$target_pos == 'same_obj'],supporto$correct_ed[supporto$target_pos == 'diff_obj'], paired = T)
cohen.d(supporto$correct_ed[supporto$target_pos == 'same_obj'],supporto$correct_ed[supporto$target_pos == 'diff_obj'], paired = T)

##### intertrial ####
####### using LAG function ########

library(dplyr)


dataset$cue_coord = (dataset$cue_pos_y == lag(dataset$cue_pos_y)) & (dataset$cue_pos_x == lag(dataset$cue_pos_x))
dataset$target_coord = (dataset$target_x == lag(dataset$target_x)) & (dataset$target_y == lag(dataset$target_y))

# Typetrial, Previous reward, target pos, prev target pos  

#selezione = data.frame(previous_target_pos = lag(dataset$target_pos), target_pos = dataset$target_pos, previous_reward = lag(dataset$reward_amount), reward = dataset$reward_amount)
selezione = data.frame(previous_target_pos = lag(dataset$target_pos), target_pos = dataset$target_pos,
                       previous_reward = lag(dataset$reward_amount), reward = dataset$reward_amount,
                       prev_outliers = lag(dataset$rt_outliers_logic), outliers = dataset$rt_outliers_logic)

selezione$keep = T
head(selezione)

# discard first trial of every block (one every 72 trials)
No_subjs = length(unique(dataset$subj_nr))
first_trial_blocks = seq(1,648*No_subjs, 72) 
selezione$keep[first_trial_blocks] = F
head(selezione) #notice first trial excluded

# outliers will be computed later
selezione$keep[selezione$reward < 0] = F # exclude wrong trials
selezione$keep[selezione$previous_reward < 0] = F # exclude trials preceded by wrong trials
selezione$keep[!selezione$outliers] = F 
selezione$keep[!selezione$prev_outliers] = F 


selezione = selezione[,c("previous_target_pos","target_pos","previous_reward", "keep")]
head(selezione)
# assign labels
selezione$previous_reward = as.character(selezione$previous_reward)
selezione[selezione == "6"] = "HIGH"
selezione[selezione == "1"] = "LOW"
selezione[selezione == "valid"] = 'V'
selezione[selezione == "same_obj"] = 'SO'
selezione[selezione == "diff_obj"] = 'DO'
selezione$typetrial = paste(selezione$previous_target_pos, "_", selezione$target_pos, sep = "") # create typetrial var
selezione$consistency = 'inconsistent'
selezione$consistency[selezione$typetrial == "V_V" | selezione$typetrial == "SO_SO" | selezione$typetrial == "DO_DO"] = "consistent"

selezione$keep[selezione$typetrial %in% c('V_V',"V_SO","V_DO","DO_V","SO_V")] = F

head(selezione)

#### add label vars to dataset ####
dataset$previous_reward = selezione$previous_reward 
dataset$previous_target_pos = selezione$previous_target_pos 
dataset$target_pos = selezione$target_pos
dataset$typetrial = selezione$typetrial
dataset$consistency = selezione$consistency
dataset$keep = selezione$keep

# we will exclude trials later after computing inv efficiency and calculating outliers
#full_dataset = dataset 
#dataset = dataset[dataset$keep,]
dataset = dataset[2:length(dataset$PROLIFIC_PID),]

rt = dataset$reaction_time
target_pos = as.factor(dataset$target_pos)
subj_nr = as.factor(dataset$number)
typetrial = as.factor(dataset$typetrial)
previous_reward = as.factor(dataset$previous_reward)
consistency = as.factor(dataset$consistency)
previous_target_pos = as.factor(dataset$previous_target_pos)
accuracy = as.logical(dataset$correct)
keep = as.logical(dataset$keep)
target_coord = as.factor(dataset$target_coord)
cue_coord = as.factor(dataset$cue_coord)

model_data = data.frame(rt, accuracy, subj_nr, previous_target_pos, target_pos, consistency, typetrial, previous_reward,target_coord,cue_coord, keep)
head(model_data)
summary(model_data)




##### violin plot ######

previous_reward
consistency
## nuovo plot 

{
  plot_data = model_data[as.logical(model_data$keep),] %>%
    group_by(subj_nr, consistency,previous_reward) %>%
    summarise(rt = mean(rt), n=n())
  
  plot_data = transform(plot_data, drwd = ifelse(previous_reward == "HIGH", 
                                                 as.numeric(consistency)-.125, # 0.25 with position_dodge = 1
                                                 as.numeric(consistency)+.125 ) )
  plot_data2 = model_data[as.logical(model_data$keep),] %>%
    group_by(subj_nr, consistency,previous_reward) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(consistency,previous_reward) %>%
    summarise(rt = mean(rt), n=n())
  
  plot_data2 = transform(plot_data2, drwd = ifelse(previous_reward == "HIGH", 
                                                   as.numeric(consistency)-.125,
                                                   as.numeric(consistency)+.125 ) )
  
  
  windows()
  p = ggplot(data=plot_data, aes(x=consistency, y=rt, fill=previous_reward))+
    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(previous_reward, consistency), shape=previous_reward),alpha= 3/10, position = position_dodge(width=0.5))+
    geom_line(aes(x=drwd, y=rt, group=interaction(subj_nr, consistency)),alpha= 3/10)+
    geom_point(data=plot_data2, aes(group=interaction(previous_reward, consistency), shape=previous_reward), size=3, position = position_dodge(width=0.5))+
    geom_line(data=plot_data2, aes(x=drwd, y=rt, group=consistency),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
} 



model_data$previous_reward
consistency
# plot RT space- and object-based effect 
{
  windows()
  p = model_data[as.logical(model_data$keep),] %>%
    group_by(subj_nr, consistency,previous_reward) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj_nr) %>%
    mutate(subj_mean = mean(rt), consistency,previous_reward) %>%
    group_by(consistency,previous_reward) %>%
    mutate(group_mean = mean(rt), subj_nr) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(consistency,previous_reward)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=consistency, y=rt, fill=previous_reward))+
    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(625,700))+
    theme_minimal()+
    #    scale_x_discrete(labels= c("different", "same"))+
    #    scale_fill_discrete(name = "Cue & Traget", labels = c("different", "same", "valid"))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(factor(dataset$subj_nr)))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Cue & Traget") + ylab("RT")
  p
}


### quick ANOVA ####

per_anova = model_data[as.logical(model_data$keep),] %>%
  group_by(subj_nr, consistency,previous_reward) %>%
  summarise(rt = mean(rt))

per_anova = ungroup(per_anova)
str(per_anova)


res.aov <- anova_test(data = per_anova, dv = rt, wid = subj_nr,  within = c(consistency, previous_reward), effect.size = "pes")
get_anova_table(res.aov)