# load packages and functions
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(viridis)
library(ggridges)
library(forcats)
library(MASS)
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'
  ) )
}                                           
plot_distr = function(...){
  # per participant.
  windows()
  p = ggplot(dataset[as.logical(dataset$correct) & dataset$rt<2400,], aes(x = rt, y = subj_nr, fill = ..x..)) +
    geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01, gradient_lwd = 1.) +
    scale_x_continuous(expand = c(0.01, 0)) +
    scale_y_discrete(expand = c(0.01, 0)) +
    scale_fill_viridis(name = "RT", option = "C") +
    labs(title = 'Reation Times per participant',
         subtitle = paste('Cut off 1700ms\n', 'Number of pp: ', length(levels(dataset$subj_nr)))) +
    theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())+
    theme_grey()
  print(p)
}
library(readr)

dataset <- read_csv("exp4_dataset.csv")

dataset = dataset[dataset$block_number!=0,]

#### assign number to participants ####
number_of_subjs = length(unique(dataset$PROLIFIC_PID)) # number of subjects in dataset
n_trials_per_subj = 72*6
dataset$subj_nr =  rep(1:number_of_subjs, each = n_trials_per_subj)
dataset$subj_nr <- as.factor(dataset$subj_nr)
length(dataset$PROLIFIC_PID)


# Rename phases
dataset$Phase = dataset$block_number == 1 | dataset$block_number == 2
dataset$Phase[dataset$Phase] = "Baseline_Bonus" 
dataset$Phase[dataset$Phase != "Baseline_Bonus" & dataset$condition == "BB"] = "Bonus"
dataset$Phase[dataset$Phase != "Baseline_Bonus" & dataset$condition == "BE"] = "Earn"
dataset$Phase = as.factor(dataset$Phase)



dataset$target_duration = as.numeric(dataset$time_delay_target) - as.numeric(dataset$time_target)
dataset$delay_duration =  as.numeric(dataset$time_mask) - as.numeric(dataset$time_delay_target)
dataset$mask_duration = as.numeric(dataset$time_keyboard_response) - as.numeric(dataset$time_mask)


# Actual reaction times
dataset$reaction_time =  as.numeric(dataset$response_time_keyboard_response) + as.numeric(dataset$delay_duration) + as.numeric(dataset$target_duration) + as.numeric(dataset$mask_duration)

summary(dataset$reaction_time)
hist(dataset$reaction_time)
# max time allowed was 1700+100+60 but doesn't change results with straight 1700ms 
#dataset = dataset[dataset$reaction_time<1700,]
# try even including everything (over 2000ms must be an error/freezing of the browser or OSWeb) doesn't change the results
dataset$reaction_time[dataset$reaction_time>2000]
dataset = dataset[dataset$reaction_time<2000,]

dataset$rt = dataset$reaction_time


# valid trials logical vector
valid_trials = dataset$target_pos == "valid"
print(length(dataset$PROLIFIC_PID)/sum(valid_trials)/number_of_subjs)

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

accuracy = dataset %>%
  group_by(subj_nr,Phase) %>%
  summarise(mean(correct))


unique(dataset$subj_nr[dataset$condition == "BB"])
unique(dataset$subj_nr[dataset$condition == "BE"])

length(unique(dataset$subj_nr[dataset$condition == "BB"]))
length(unique(dataset$subj_nr[dataset$condition == "BE"]))


plot_pp_outliers() # 11,38,5,19,3
dataset = dataset[!dataset$subj_nr %in% c(11,38,5,19,3),]
length(unique(dataset$subj_nr[dataset$condition == "BE"]))
length(unique(dataset$subj_nr[dataset$condition == "BB"]))

# Need to reject another participant for balancing: apply again the plot_pp_outliers() and reject the first "outlier" participant in the BB group
plot_pp_outliers() # 8, 20, 50. 8 is in the BB group -> reject participant 8.
dataset = dataset[!dataset$subj_nr %in% c(8),]

dataset$subj_nr = factor(dataset$subj_nr)


length(levels(dataset$subj_nr))
acc = dataset %>%
  group_by(subj_nr)%>%
  summarise(accuracy = mean(correct))
summary(acc$accuracy)



#### create rew trg rel var ####
# object based effect and reward relation
{
  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)
  obj_eff_rwd_data$subj = factor(obj_eff_rwd_data$subj)
  
}

# look at n
obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
  group_by(rwd, trg,subj, Phase) %>%
  summarise(rt = mean(rt), n=n()) %>%
  group_by(rwd, trg, Phase) %>%
  summarise(rt = mean(rt), n=round(mean(n))) %>% print(n = 16)


# 
# dataset[dataset$condition == "BB",] %>%
#   group_by(PROLIFIC_PID) %>%
#   summarise((sum(reward_amount)*0.001))
# 
# dataset[dataset$condition == "BE",] %>%
#   subset(block_number == "1" | block_number == "2") %>%             
#   group_by(PROLIFIC_PID) %>%
#   summarise((sum(reward_amount)*0.001)*3)




#### plot obj-based eff and reward ####

# put together BB & BE for the first two blocks

baseline_bonus = obj_eff_rwd_data[obj_eff_rwd_data$block_number == 1 | obj_eff_rwd_data$block_number == 2,]
{
  windows()
  p = baseline_bonus[as.logical(baseline_bonus$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, trg, n) %>%
    group_by(rwd, trg) %>%
    mutate(group_mean = mean(rt), subj, n) %>%
    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,940))+
    scale_y_continuous(breaks=seq(600, 940, 20))+
    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("First Learning Phase (block 1 & 2 all participants)\n(error bar: intra-subj variability)\nNumber of participants per between condition",length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"]))))+ 
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}

# plot BB & BE for blocks after two

BB_BE = obj_eff_rwd_data[obj_eff_rwd_data$block_number > 2,]



# BB & BE
{
  windows()
  p = BB_BE[as.logical(BB_BE$correct),] %>%
    group_by(rwd, trg,subj, condition) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(subj, condition) %>%
    mutate(subj_mean = mean(rt), rwd, trg, n) %>%
    group_by(rwd, trg, condition) %>%
    mutate(group_mean = mean(rt), subj, n) %>%
    group_by(condition) %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, trg,  condition)%>%
    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(570,800))+
    scale_y_continuous(breaks=seq(570, 800, 20))+
    theme_minimal()+
    facet_grid(. ~ condition)+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+ 
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    ggtitle(paste("BB vs BE\n(error bar: intra-subj variability)\nNumber of participants per between condition",length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"]))))+ 
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}


# Baseline bonus + BB_BE
{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,subj, Phase, condition) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(subj, Phase, condition) %>%
    mutate(subj_mean = mean(rt), rwd, trg, n) %>%
    group_by(rwd, trg, Phase,condition) %>%
    mutate(group_mean = mean(rt), subj, n) %>%
    group_by(condition) %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, trg, Phase, condition)%>%
    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,940))+
    scale_y_continuous(breaks=seq(600, 940, 20))+
    theme_minimal()+
    facet_grid(condition ~ Phase)+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+ 
    scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
    ggtitle(paste("BB vs BE\n(error bar: intra-subj variability)\nNumber of participants per between condition",length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"]))))+ 
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}

unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"])
unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BE"])

length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"]))
length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BE"]))


{
  windows()
  p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>%
    group_by(rwd, trg,condition,Phase,subj) %>%
    summarise(rt = mean(rt)) %>%             
    group_by(rwd, trg,condition,Phase) %>%
    summarise(rt = mean(rt)) %>%
    ggplot(aes(x=rwd, y=rt, fill=trg))+
    geom_bar(stat="identity", width=0.5, position = "dodge")+
    coord_cartesian(ylim = c(580,900))+
    theme_minimal()+
    #    scale_x_discrete(labels= c("Low reward", "Cue not rew ass", "Cue rew ass"))+
    #    scale_fill_discrete(name = "Type trial", 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("Trial type") + ylab("RT")+
    facet_grid(Phase ~ condition)
  p
}

######## STATS INVALID TRIALS #########
#### 2X2X2 trg X rwd - 2nd Phase - no low rewrd ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase != "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# without low reward baseline

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

###ACC####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg, Phase)%>%
  summarise(acc = mean(correct), n=n())

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


windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(acc), rwd, trg, Phase, n) %>%
  group_by(rwd,Phase, trg) %>%
  mutate(group_mean = mean(acc), subj, n) %>%
  group_by() %>%
  mutate(isv = (acc - subj_mean + group_mean)) %>%
  group_by(rwd,Phase, trg)%>%
  summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  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(0.7,1))+
  scale_y_continuous(breaks=seq(0.7,1, 0.05))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Accuracy")+
  geom_text(aes(label=round(acc,3)), position=position_dodge(width=0.9), vjust=-0.25)+
  facet_grid(. ~Phase)


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

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



per_anova_earn %>%
  group_by(rwd,trg, Phase)%>%
  summarise(acc = mean(acc)) %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")

# main effect of reward
per_anova_earn %>%
  group_by(trg)%>%
  summarise(acc = mean(acc))

# effect of reward X phase
per_anova_earn %>%
  group_by(Phase)%>%
  summarise(acc = mean(acc))

library(lme4)
library(lmerTest)
options(contrasts = c("contr.sum","contr.poly")) 

ciaobelli = glmer(correct ~ rwd * trg * Phase + (1 | subj),
                  data = obj_eff_2ndPhase_hr,
                  family = binomial(link = "logit"))     #The link-function describing the relationship between your  linear model and the binomial distribution of your outcome

r_int<- ranef(ciaobelli)$subj$`(Intercept)`
qqnorm(r_int)
qqline(r_int)
shapiro.test(r_int)

summary(ciaobelli)


### RT ####
## Statistical Model

per_anova_earn = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
  group_by(subj, rwd,trg, Phase)%>%
  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,
  between = Phase, within = c(rwd, trg),
  effect.size = "pes")

get_anova_table(res.aov)
set.seed(1000)
bf = anovaBF(rt ~ rwd*trg*Phase + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 100000)
#save(bf, file = "BF_2X2X2_trgXrwd_2ndPhase_nolowrewrd.RData")

against_best = bf/bf

per_anova_earn %>%
  group_by(Phase,rwd,trg)%>%
  summarise(rt = mean(rt)) 

#### inv eff ####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg, Phase)%>%
  summarise(acc = mean(correct), n = n())

supporto = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
  group_by(subj, rwd,trg, Phase)%>%
  summarise(ciao = mean(rt), n=n())

per_anova_earn$ciao = supporto$ciao
per_anova_earn$inv_eff = per_anova_earn$ciao / per_anova_earn$acc


per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)
# hist(per_anova_earn$inv_eff)
# hist(per_anova_earn$acc)
# hist(per_anova_earn$rt)


res.aov <- anova_test(data = per_anova_earn, 
                      dv = inv_eff, 
                      wid = subj,
                      between = "Phase",
                      within = c(rwd, trg), 
                      effect.size = "pes")
get_anova_table(res.aov)
# equivalence test #####
calc_LL_LE = function(z){
  LL = (
    z$inv_eff[
      z$Phase == "Bonus" &
        z$rwd == "congruent" &
        z$trg == "diff_obj" 
    ] -
      z$inv_eff[
        z$Phase == "Bonus" &
          z$rwd == "congruent" &
          z$trg == "same_obj" 
      ] 
  )-
    (
      z$inv_eff[
        z$Phase == "Bonus" &
          z$rwd == "incongruent" &
          z$trg == "diff_obj" 
      ] -
        z$inv_eff[
          z$Phase == "Bonus" &
            z$rwd == "incongruent" &
            z$trg == "same_obj" 
        ] 
    )
  
  LE = (
    z$inv_eff[
      z$Phase == "Earn" &
        z$rwd == "congruent" &
        z$trg == "diff_obj" 
    ] -
      z$inv_eff[
        z$Phase == "Earn" &
          z$rwd == "congruent" &
          z$trg == "same_obj" 
      ] 
  )-
    (
      z$inv_eff[
        z$Phase == "Earn" &
          z$rwd == "incongruent" &
          z$trg == "diff_obj" 
      ] -
        z$inv_eff[
          z$Phase == "Earn" &
            z$rwd == "incongruent" &
            z$trg == "same_obj" 
        ] 
    )
  
  return(list(LL,LE))
}


LL_LE = calc_LL_LE(per_anova_earn)

library(simpleboot)
LL <- data.frame(group = 'LL', inv_eff = LL_LE[[1]])
LE <- data.frame(group = 'LE', inv_eff = LL_LE[[2]])
samples <- rbind(LL,LE)
ggplot(samples,aes(x = inv_eff, col = group, fill = group, group = group)) + geom_density(alpha=0.2)


bootstrapped <- two.boot(LL$inv_eff, LE$inv_eff, mean, 100000)

bootstrapped_mean_diff <- data.frame(bootstrapped$t)
colnames(bootstrapped_mean_diff) <- 'mean_diffs'

ggplot(bootstrapped_mean_diff, aes(x=mean_diffs)) + 
  geom_histogram(bins=200, alpha=0.8) +
  geom_vline(xintercept = bootstrapped$t0, size = 1.5) +
  geom_vline(xintercept = quantile(bootstrapped_mean_diff$mean_diffs, 0.05), size = 1, linetype = 'dashed') +
  geom_vline(xintercept = quantile(bootstrapped_mean_diff$mean_diffs, 0.95), size = 1, linetype = 'dashed') +
  xlab('sample mean') + ylab('bootstrapped count')

equivalence_bound = 15/.75 # 15ms / "average" performance to get IE score 
t.test(LL$inv_eff-equivalence_bound,LE$RT,alternative = "greater")
t.test(LL$inv_eff+equivalence_bound,LE$RT,alternative = "greater")


library("TOSTER")
TOSTtwo(m1=100.64,m2=100.48,sd1=14.1,sd2=14.9,n1=39343,n2=40033,low_eqbound_d=-0.05, high_eqbound_d=0.05, alpha = 0.05, var.equal=FALSE)

#### t test - 1st phase - low reward ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# low reward baseline

obj_eff_2ndPhase_hr = obj_eff_2ndPhase[obj_eff_2ndPhase$side_high_reward == "to_be_assigned",]
obj_eff_2ndPhase_hr$rwd = factor(obj_eff_2ndPhase_hr$rwd)  # relevel with only high rwd trials


## Statistical Model

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

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

t_test(per_anova_earn, rt ~ trg, paired = T)


per_anova_earn %>%
  group_by(trg)%>%
  summarise(mean(rt))

library(effsize)
cohens_d(per_anova_earn, rt ~ trg, paired = T)






#### 2X2 trg X rwd - 1st Phase - no low rewrd ####
obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# without low reward baseline

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


###ACC####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct), n=n())

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

windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(acc), rwd, trg, n) %>%
  group_by(rwd, trg) %>%
  mutate(group_mean = mean(acc), subj, n) %>%
  group_by() %>%
  mutate(isv = (acc - subj_mean + group_mean)) %>%
  group_by(rwd, trg)%>%
  summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  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(0.7,1))+
  scale_y_continuous(breaks=seq(0.7,1, 0.05))+
  scale_x_discrete(labels= c("Target on high reward", "Traget on low reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Accuracy")+
  geom_text(aes(label=round(acc,3)), position=position_dodge(width=0.9), vjust=-0.25)



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

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

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

# 
# 
# library(lme4)
# library(lmerTest)
# options(contrasts = c("contr.sum","contr.poly")) 
# 
# ciaobelli = glmer(correct ~ rwd * trg + (1 | subj),
#                   data = obj_eff_2ndPhase_hr,
#                   family = binomial(link = "logit"))     #The link-function describing the relationship between your  linear model and the binomial distribution of your outcome
# 
# r_int<- ranef(ciaobelli)$subj$`(Intercept)`
# qqnorm(r_int)
# qqline(r_int)
# shapiro.test(r_int)
# 
# summary(ciaobelli)


#### inv eff ####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct), n = n())

supporto = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(ciao = mean(rt), n=n())

per_anova_earn$ciao = supporto$ciao
per_anova_earn$inv_eff = per_anova_earn$ciao / per_anova_earn$acc


per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)
# hist(per_anova_earn$inv_eff)
# hist(per_anova_earn$acc)
# hist(per_anova_earn$rt)


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)

per_anova_earn %>%
  group_by(rwd)%>%
  summarise(inv_eff = mean(inv_eff), n=n())

per_anova_earn %>%
  group_by(trg)%>%
  summarise(inv_eff = mean(inv_eff), n=n())

windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(inv_eff), rwd, trg, n) %>%
  group_by(rwd, trg) %>%
  mutate(group_mean = mean(inv_eff), subj, n) %>%
  group_by() %>%
  mutate(isv = (inv_eff - subj_mean + group_mean)) %>%
  group_by(rwd, trg)%>%
  summarise(inv_eff = mean(inv_eff), isv = sem_fun(isv))  %>%
  ggplot(aes(y=inv_eff,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")+
  geom_errorbar(aes(ymin=inv_eff-isv, ymax=inv_eff+isv), width=.2,position=position_dodge(.5))+ 
  coord_cartesian(ylim = c(600,1200))+
  scale_y_continuous(breaks=seq(600,1200, 100))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Inverse efficiency score")+
  geom_text(aes(label=round(inv_eff)), position=position_dodge(width=0.9), vjust=-0.25)



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

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


windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(inv_eff_ed), rwd, trg, n) %>%
  group_by(rwd, trg) %>%
  mutate(group_mean = mean(inv_eff_ed), subj, n) %>%
  group_by() %>%
  mutate(isv = (inv_eff_ed - subj_mean + group_mean)) %>%
  group_by(rwd, trg)%>%
  summarise(inv_eff_ed = mean(inv_eff_ed), isv = sem_fun(isv))  %>%
  ggplot(aes(y=inv_eff_ed,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")+
  geom_errorbar(aes(ymin=inv_eff_ed-isv, ymax=inv_eff_ed+isv), width=.2,position=position_dodge(.5))+ 
  coord_cartesian(ylim = c(1.260,1.266))+
  #  scale_y_continuous(breaks=seq(0.69706,0.69719, 0.00005))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Inverse efficiency score - Tranformed values")



###RT####
## Statistical Model

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

#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)
{
  windows()
  p = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
    group_by(rwd, trg,subj) %>%
    summarise(rt = mean(rt), n=n()) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, trg, n) %>%
    group_by(rwd, trg) %>%
    mutate(group_mean = mean(rt), subj, n) %>%
    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,940))+
    scale_y_continuous(breaks=seq(600, 940, 20))+
    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("First Learning Phase (block 1 & 2 all participants)\n(error bar: intra-subj variability)\nNumber of participants per between condition",length(unique(obj_eff_rwd_data$subj[obj_eff_rwd_data$condition == "BB"]))))+ 
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}

set.seed(1000)
bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 100000)

#save(bf, file = "BF_2X2_trgXrwd_1stPhase_nolowrewrd.RData")

options(contrasts = c("contr.sum","contr.poly")) 
model = lmer(rt ~ rwd*trg + (1|subj),
             data =  obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),])
summary(model)
anova(model)
windows()
hist(resid(model))
windows()
qqnorm(resid(model))
qqline(resid(model))

# Gamma(link = "identity") leads to better looking residuals, results are the same
model = glmer(rt ~ rwd*trg + (1|subj), data = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),], family = Gamma(link = "identity"))
summary(model)
windows()
hist(resid(model))
windows()
qqnorm(resid(model))
qqline(resid(model))

# look at n
obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
  group_by(rwd, trg,subj) %>%
  summarise(rt = mean(rt), n=n()) %>%
  group_by(rwd, trg) %>%
  summarise(rt = mean(rt), n=round(mean(n))) %>% print(n = 16)

#### add group variable ####
# add group variable to see if it makes any difference, that is, if the groups are different to start with

###ACC####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg, condition)%>%
  summarise(acc = mean(correct), n=n())

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

windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(acc), rwd, trg, condition, n) %>%
  group_by(rwd, trg, condition) %>%
  mutate(group_mean = mean(acc), subj, n) %>%
  group_by() %>%
  mutate(isv = (acc - subj_mean + group_mean)) %>%
  group_by(rwd, trg, condition)%>%
  summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  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(0.7,1))+
  scale_y_continuous(breaks=seq(0.7,1, 0.05))+
  scale_x_discrete(labels= c("Target on high reward", "Traget on low reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Accuracy")+
  geom_text(aes(label=round(acc,3)), position=position_dodge(width=0.9), vjust=-0.25)+
  facet_grid(. ~ condition)



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

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


windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(acc_ed), rwd, trg, condition, n) %>%
  group_by(rwd, trg, condition) %>%
  mutate(group_mean = mean(acc_ed), subj, n) %>%
  group_by() %>%
  mutate(isv = (acc_ed - subj_mean + group_mean)) %>%
  group_by(rwd, trg, condition)%>%
  summarise(acc_ed = mean(acc_ed), isv = sem_fun(isv))  %>%
  ggplot(aes(y=acc_ed,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")+
  geom_errorbar(aes(ymin=acc_ed-isv, ymax=acc_ed+isv), width=.2,position=position_dodge(.5))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Accuracy - Tranformed")+
  facet_grid(. ~ condition)




#### inv eff ####
per_anova_earn = obj_eff_2ndPhase_hr %>%
  group_by(subj, rwd,trg,condition)%>%
  summarise(acc = mean(correct), n = n())

supporto = obj_eff_2ndPhase_hr[as.logical(obj_eff_2ndPhase_hr$correct),] %>%
  group_by(subj, rwd,trg,condition)%>%
  summarise(ciao = mean(rt), n=n())

per_anova_earn$ciao = supporto$ciao
per_anova_earn$inv_eff = per_anova_earn$ciao / per_anova_earn$acc


per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)
# hist(per_anova_earn$inv_eff)
# hist(per_anova_earn$acc)
# hist(per_anova_earn$rt)


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

windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(inv_eff), rwd, trg, n) %>%
  group_by(rwd, trg) %>%
  mutate(group_mean = mean(inv_eff), subj, n) %>%
  group_by() %>%
  mutate(isv = (inv_eff - subj_mean + group_mean)) %>%
  group_by(rwd, trg)%>%
  summarise(inv_eff = mean(inv_eff), isv = sem_fun(isv))  %>%
  ggplot(aes(y=inv_eff,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")+
  geom_errorbar(aes(ymin=inv_eff-isv, ymax=inv_eff+isv), width=.2,position=position_dodge(.5))+ 
  coord_cartesian(ylim = c(600,1200))+
  scale_y_continuous(breaks=seq(600,1200, 100))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Inverse efficiency score")+
  geom_text(aes(label=round(inv_eff)), position=position_dodge(width=0.9), vjust=-0.25)



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

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


windows()
per_anova_earn %>%
  group_by(subj) %>%
  mutate(subj_mean = mean(inv_eff_ed), rwd, trg, n) %>%
  group_by(rwd, trg) %>%
  mutate(group_mean = mean(inv_eff_ed), subj, n) %>%
  group_by() %>%
  mutate(isv = (inv_eff_ed - subj_mean + group_mean)) %>%
  group_by(rwd, trg)%>%
  summarise(inv_eff_ed = mean(inv_eff_ed), isv = sem_fun(isv))  %>%
  ggplot(aes(y=inv_eff_ed,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")+
  geom_errorbar(aes(ymin=inv_eff_ed-isv, ymax=inv_eff_ed+isv), width=.2,position=position_dodge(.5))+ 
  coord_cartesian(ylim = c(1.260,1.266))+
  #  scale_y_continuous(breaks=seq(0.69706,0.69719, 0.00005))+
  scale_x_discrete(labels= c("Low reward", "High reward"))+ 
  scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+
  ggtitle("Inverse efficiency score - Tranformed values")



###RT####
## Statistical Model

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

#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),
  between = "condition",
  effect.size = "pes",
  detailed = T)

get_anova_table(res.aov)




#### t test - 1st phase - low reward ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# low reward baseline

obj_eff_2ndPhase_hr = obj_eff_2ndPhase[obj_eff_2ndPhase$side_high_reward == "to_be_assigned",]
obj_eff_2ndPhase_hr$rwd = factor(obj_eff_2ndPhase_hr$rwd)  # relevel with only high rwd trials


## Statistical Model

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

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

t_test(per_anova_earn, rt ~ trg, paired = T)


per_anova_earn %>%
  group_by(trg)%>%
  summarise(mean(rt))

library(effsize)
cohens_d(per_anova_earn, rt ~ trg, paired = T)






#### t test - 1st phase - incongruent ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# incongruent only

obj_eff_incong = obj_eff_2ndPhase[obj_eff_2ndPhase$rwd == "incongruent",]
obj_eff_incong$rwd = factor(obj_eff_incong$rwd)  # relevel with only high rwd trials


## Statistical Model

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

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

t_test(per_anova_earn, rt ~ trg, paired = T)


per_anova_earn %>%
  group_by(trg)%>%
  summarise(mean(rt))

library(effsize)
cohens_d(per_anova_earn, rt ~ trg, paired = T)



#### 2X2 trg X rwd - 2nd Phase - only Bonus phase - no low rewrd ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase != "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# without low reward baseline

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


obj_eff_Bonus = obj_eff_2ndPhase_hr[obj_eff_2ndPhase_hr$condition == "BB",]
obj_eff_Bonus$condition = factor(obj_eff_Bonus$condition) 

###ACC####
per_anova_earn = obj_eff_Bonus %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct))

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

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

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



per_anova_earn %>%
  group_by(rwd,trg)%>%
  summarise(acc = mean(acc)) %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")



library(lme4)
library(lmerTest)
options(contrasts = c("contr.sum","contr.poly")) 

ciaobelli = glmer(correct ~ rwd * trg + (1 | subj),
                  data = obj_eff_Bonus,
                  family = binomial(link = "logit"))     #The link-function describing the relationship between your  linear model and the binomial distribution of your outcome

r_int<- ranef(ciaobelli)$subj$`(Intercept)`
qqnorm(r_int)
qqline(r_int)
shapiro.test(r_int)

summary(ciaobelli)


####RT####
## Statistical Model

per_anova_earn = obj_eff_Bonus[as.logical(obj_eff_Bonus$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)

per_anova_earn %>%
  group_by(rwd)%>%
  summarise(rt = mean(rt))


per_anova_earn %>%
  group_by(trg)%>%
  summarise(rt = mean(rt))



set.seed(1000)
bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 100000)

against_best = bf/bf
against_best[3,4]

options(contrasts = c("contr.sum","contr.poly")) 
model = lmer(rt ~ rwd*trg + (1|subj),
             data = obj_eff_Bonus[as.logical(obj_eff_Bonus$correct),])
summary(model)
anova(model)
windows()
hist(resid(model))
windows()
qqnorm(resid(model))
qqline(resid(model))

# Gamma(link = "identity") leads to better looking residuals, results are the same
model = glmer(rt ~ rwd*trg + (1|subj), data = obj_eff_Bonus[as.logical(obj_eff_Bonus$correct),], family = Gamma(link = "identity"))
summary(model)
windows()
hist(resid(model))
windows()
qqnorm(resid(model))
qqline(resid(model))

# look at n
obj_eff_Bonus[as.logical(obj_eff_Bonus$correct),] %>%
  group_by(rwd, trg,subj) %>%
  summarise(rt = mean(rt), n=n()) %>%
  group_by(rwd, trg) %>%
  summarise(rt = mean(rt), n=round(mean(n))) %>% print(n = 16)

#### inv eff ####
per_anova_earn = obj_eff_Bonus %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct), n = n())

supporto = obj_eff_Bonus[as.logical(obj_eff_Bonus$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(ciao = mean(rt), n=n())

per_anova_earn$ciao = supporto$ciao
per_anova_earn$inv_eff = per_anova_earn$ciao / per_anova_earn$acc


per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)
# hist(per_anova_earn$inv_eff)
# hist(per_anova_earn$acc)
# hist(per_anova_earn$rt)


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 trg X rwd - 2nd Phase - only Extinction phase - no low rewrd ####


obj_eff_2ndPhase = obj_eff_rwd_data[obj_eff_rwd_data$Phase != "Baseline_Bonus",]
obj_eff_2ndPhase$Phase = factor(obj_eff_2ndPhase$Phase)  # relevel


# without low reward baseline

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

# only Bonus phase

obj_eff_Bonus = obj_eff_2ndPhase_hr[obj_eff_2ndPhase_hr$condition == "BE",]
obj_eff_Bonus$condition = factor(obj_eff_Bonus$condition)  # relevel with only high rwd trials

###ACC####
per_anova_earn = obj_eff_Bonus %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct))

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

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

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

per_anova_earn %>%
  group_by(rwd)%>%
  summarise(acc = mean(acc))

per_anova_earn %>%
  group_by(trg)%>%
  summarise(acc = mean(acc))


per_anova_earn %>%
  group_by(rwd,trg)%>%
  summarise(acc = mean(acc)) %>%
  ggplot(aes(y=acc,x=rwd,fill=trg))+
  geom_bar(stat="identity", width=0.5, position = "dodge")



library(lme4)
library(lmerTest)
options(contrasts = c("contr.sum","contr.poly")) 

ciaobelli = glmer(correct ~ rwd * trg + (1 | subj),
                  data = obj_eff_Bonus,
                  family = binomial(link = "logit"))     #The link-function describing the relationship between your  linear model and the binomial distribution of your outcome

r_int<- ranef(ciaobelli)$subj$`(Intercept)`
qqnorm(r_int)
qqline(r_int)
shapiro.test(r_int)

summary(ciaobelli)


####RT####
## Statistical Model

per_anova_earn = obj_eff_Bonus[as.logical(obj_eff_Bonus$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)

per_anova_earn %>%
  group_by(rwd)%>%
  summarise(rt = mean(rt))

per_anova_earn %>%
  group_by(trg)%>%
  summarise(rt = mean(rt))


set.seed(1000)
bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, 
             whichRandom="subj", iterations = 100000)

against_best = bf/bf
against_best[3,4]


#### inv eff ####
per_anova_earn = obj_eff_Bonus %>%
  group_by(subj, rwd,trg)%>%
  summarise(acc = mean(correct), n = n())

supporto = obj_eff_Bonus[as.logical(obj_eff_Bonus$correct),] %>%
  group_by(subj, rwd,trg)%>%
  summarise(ciao = mean(rt), n=n())

per_anova_earn$ciao = supporto$ciao
per_anova_earn$inv_eff = per_anova_earn$ciao / per_anova_earn$acc


per_anova_earn = ungroup(per_anova_earn)
str(per_anova_earn)
# hist(per_anova_earn$inv_eff)
# hist(per_anova_earn$acc)
# hist(per_anova_earn$rt)


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)

######################### VALID #############################

# condition becomes group
dataset$group = factor(dataset$condition)

# Phase becomes phase
dataset$phase = dataset$Phase %in% c("Bonus", "Earn")
dataset$phase = as.numeric(dataset$phase)
dataset$phase = factor(dataset$phase)

valid_rwd_data = dataset[dataset$target_pos == "valid",] # taking into account only invalid trials

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

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

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

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

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

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

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


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




acc = dataset %>%
  group_by(PROLIFIC_PID)%>%
  summarise(accuracy = mean(correct))
summary(acc$accuracy)
# number of trials in the smallest cell per subj
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  group_by(rwd, trg,subj,Phase) %>%
  summarise(rt = mean(rt), n=n()) %>%
  group_by(rwd, trg, Phase) %>%
  summarise(rt = mean(rt), n=mean(n)) %>% print(n = 100)

# group stats ####
# in general no conditional difference
dataset[as.logical(dataset$correct),] %>%
  subset(Phase == "Baseline_Bonus") %>%
  group_by(subj_nr, condition)%>%
  summarise(rt = mean(rt))%>%
  group_by(condition)%>%
  summarise( sd = sd(rt), rt = mean(rt))

# validly cue trials
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "0") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  group_by(group)%>%
  summarise( sd = sd(rt), rt = mean(rt))

# plot RT, rew-trg and cue-trg 
{
  windows()
  p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
    group_by(rwd, subj, phase, group) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, phase,group) %>%
    group_by(rwd,  phase,group) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd,  phase,group)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=rt, fill=rwd))+
    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(550,880))+
    theme_minimal()+
    scale_x_discrete(labels= c("Baseline", "High reward", "Low reward"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Valid"))+
    scale_y_continuous(breaks=seq(550, 880, 50))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(valid_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")+
    facet_grid(group ~ phase)
  p
  }



############### STATS VALID TRIALS ###############
# phase X group X rwd ######

## RT ####

per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  group_by(subj, rwd, phase, group)%>%
  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,
  between = group, within = c(rwd, phase),
  effect.size = "pes")

get_anova_table(res.aov)




# 1st phase:  rwd  #####
{
  windows()
  p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
    subset(phase == "0") %>%
    group_by(rwd, subj) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd) %>%
    group_by(rwd) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=rt, fill=rwd))+
    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(550,880))+
    theme_minimal()+
    scale_x_discrete(labels= c("Baseline", "High reward", "Low reward"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Valid"))+
    scale_y_continuous(breaks=seq(550, 880, 50))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(valid_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}



{
  windows()
  p = valid_rwd_data %>%
    subset(phase == "0") %>%
    group_by(rwd,subj) %>%
    summarise(acc = mean(correct)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(acc), rwd ) %>%
    group_by(rwd) %>%
    mutate(group_mean = mean(acc), subj) %>%
    group_by() %>%
    mutate(isv = (acc - subj_mean + group_mean)) %>%
    group_by(rwd)%>%
    summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=acc, fill=rwd))+
    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(0.75, 1))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Cue & Target", labels = c("Valid"))+
    scale_y_continuous(breaks=seq(0.75, 1, 0.05))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(valid_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("ACC")
  
  p
}

###RT####
## Statistical Model

per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "0") %>%
  group_by(subj, rwd)%>%
  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),
  effect.size = "pes")

get_anova_table(res.aov)


pwc <- per_anova_earn %>%
  pairwise_t_test(
    rt ~ rwd, paired = TRUE,
    p.adjust.method = "bonferroni",
    detailed = T
  )
data.frame(pwc)

library(effsize)
cohen.d(per_anova_earn$rt[per_anova_earn$rwd == "low_reward"],per_anova_earn$rt[per_anova_earn$rwd == "congruent"], paired = T)
cohen.d(per_anova_earn$rt[per_anova_earn$rwd == "incongruent"],per_anova_earn$rt[per_anova_earn$rwd == "congruent"], paired = T)
cohen.d(per_anova_earn$rt[per_anova_earn$rwd == "incongruent"],per_anova_earn$rt[per_anova_earn$rwd == "low_reward"], paired = T)


valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "0") %>%
  group_by(subj, rwd)%>%
  summarise(rt = mean(rt))%>%
  group_by(rwd)%>%
  summarise(rt = mean(rt))

### ACC ####

per_anova_earn = valid_rwd_data %>%
  subset(phase == "0")%>%
  group_by(subj,rwd)%>%
  summarise(acc = mean(correct))

#str(per_anova_earn)
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),
                      effect.size = "pes")
get_anova_table(res.aov, correction = "GG")


pwc <- per_anova_earn %>%
  pairwise_t_test(
    correct_ed ~ rwd, paired = TRUE,
    p.adjust.method = "bonferroni",
    detailed = T
  )
data.frame(pwc)

library(effsize)
cohen.d(per_anova_earn$correct_ed[per_anova_earn$rwd == "low_reward"],per_anova_earn$correct_ed[per_anova_earn$rwd == "congruent"], paired = T)
cohen.d(per_anova_earn$correct_ed[per_anova_earn$rwd == "incongruent"],per_anova_earn$correct_ed[per_anova_earn$rwd == "congruent"], paired = T)
cohen.d(per_anova_earn$correct_ed[per_anova_earn$rwd == "incongruent"],per_anova_earn$correct_ed[per_anova_earn$rwd == "low_reward"], paired = T)



valid_rwd_data %>%
  subset(phase == "0")%>%
  group_by(subj,rwd)%>%
  summarise(acc = mean(correct))%>%
  group_by(rwd)%>%
  summarise(acc = mean(acc))




# 2nd phase, both:  rwd  #####
{
  windows()
  p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
    subset(phase == "1") %>%
    group_by(rwd, subj, group) %>%
    summarise(rt = mean(rt)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(rt), rwd, group) %>%
    group_by(rwd, group) %>%
    mutate(group_mean = mean(rt), subj) %>%
    group_by() %>%
    mutate(isv = (rt - subj_mean + group_mean)) %>%
    group_by(rwd, group)%>%
    summarise(rt = mean(rt), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=rt, fill=group))+
    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(550,880))+
    theme_minimal()+
    scale_x_discrete(labels= c("Baseline", "High reward", "Low reward"))+
    scale_fill_discrete(name = "Group", labels = c("LL", "LE"))+
    scale_y_continuous(breaks=seq(550, 880, 50))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(valid_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("RT")
  p
}



{
  windows()
  p = valid_rwd_data %>%
    subset(phase == "1") %>%
    group_by(rwd, subj, group) %>%
    summarise(acc = mean(correct)) %>%
    group_by(subj) %>%
    mutate(subj_mean = mean(acc), rwd, group ) %>%
    group_by(rwd, group) %>%
    mutate(group_mean = mean(acc), subj) %>%
    group_by() %>%
    mutate(isv = (acc - subj_mean + group_mean)) %>%
    group_by(rwd, group)%>%
    summarise(acc = mean(acc), isv = sem_fun(isv))  %>%
    ggplot(aes(x=rwd, y=acc, fill=group))+
    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(0.75, 1))+
    theme_minimal()+
    scale_x_discrete(labels= c("Low reward", "congruent", "incongruent"))+
    scale_fill_discrete(name = "Group", labels = c("LL", "LE"))+
    scale_y_continuous(breaks=seq(0.75, 1, 0.05))+
    ggtitle(paste("Object-based effect\n(error bar: intra-subj variability)\nNumber of participants",length(levels(valid_rwd_data$subj))))+
    theme(plot.title = element_text(hjust = 0.5))+
    xlab("Reward & Target") + ylab("ACC")
  
  p
}

###RT####
## Statistical Model

per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1") %>%
  group_by(subj, rwd, group)%>%
  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),
  between = group,
  effect.size = "pes")

get_anova_table(res.aov)

valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1") %>%
  group_by(subj, rwd, group)%>%
  summarise(rt = mean(rt))%>%
  group_by(group)%>%
  summarise(rt = mean(rt))

valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" ) %>%
  group_by(subj, group, rwd)%>%
  summarise(rt = mean(rt))%>%
  group_by(rwd,group)%>%
  summarise(rt = mean(rt))

library(rstatix)
# LOW REWARD - LL vs LE
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "incongruent") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  t_test(rt ~ group)
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "incongruent") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  cohens_d(rt ~ group)

# BASELINE - LL vs LE
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "low_reward") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  t_test(rt ~ group)
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "low_reward") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  cohens_d(rt ~ group)

# HIGH REWARD - LL vs LE
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "congruent") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  t_test(rt ~ group)
valid_rwd_data[as.logical(valid_rwd_data$correct),] %>%
  subset(phase == "1" & rwd == "congruent") %>%
  group_by(subj, group)%>%
  summarise(rt = mean(rt))%>%
  ungroup()%>%
  cohens_d(rt ~ group)


### ACC ####

per_anova_earn = valid_rwd_data %>%
  subset(phase == "1")%>%
  group_by(subj, rwd, group)%>%
  summarise(acc = mean(correct))

#str(per_anova_earn)
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),
                      between = group,
                      effect.size = "pes")
get_anova_table(res.aov, correction = "GG")



valid_rwd_data %>%
  subset(phase == "1")%>%
  group_by(subj,rwd, group)%>%
  summarise(acc = mean(correct))%>%
  group_by(group)%>%
  summarise(acc = mean(acc))

valid_rwd_data %>%
  subset(phase == "1")%>%
  group_by(subj,rwd, group)%>%
  summarise(acc = mean(correct))%>%
  group_by(rwd,group)%>%
  summarise(acc = mean(acc))

library(rstatix)
# LOW REWARD - LL vs LE
valid_rwd_data%>%
  subset(phase == "1" & rwd == "incongruent") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  t_test(correct ~ group)
valid_rwd_data%>%
  subset(phase == "1" & rwd == "incongruent") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  cohens_d(correct ~ group)

# BASELINE - LL vs LE
valid_rwd_data%>%
  subset(phase == "1" & rwd == "low_reward") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  t_test(correct ~ group)
valid_rwd_data%>%
  subset(phase == "1" & rwd == "low_reward") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  cohens_d(correct ~ group)

# HIGH REWARD - LL vs LE
valid_rwd_data%>%
  subset(phase == "1" & rwd == "congruent") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  t_test(correct ~ group)
valid_rwd_data%>%
  subset(phase == "1" & rwd == "congruent") %>%
  group_by(subj, group)%>%
  summarise(correct = mean(correct))%>%
  ungroup()%>%
  cohens_d(correct ~ group)