# 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("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) } options(contrasts = c("contr.sum","contr.poly")) prepare_data = function(path_to_dataset = "~/OneDrive - University of Birmingham/Desktop/Motivated Cognition Lab/Coloured rectangles + reward 2/Results/Pilot_subj.csv", suppressPlot = FALSE){ # load data, show some general stats and graphs #### pakages #### # source("IQR_custom.R") # function for Inter-Quartile Range - not needed for now, we have 1700ms answer limit library(dplyr) library(ggplot2) sem_fun = function(x) sd(x)/sqrt(length(x)) #sem library(readr) #### Load dataset #### dataset <- read_csv(path_to_dataset) if (!suppressPlot){ View(dataset) } number_of_trials = 72*10 dataset = as.data.frame(dataset) # number of subjects in dataset number_of_subjs = length(unique(dataset$PROLIFIC_PID)) # delete practice and print summary and boxplot of number of times participants took practice phase dataset = dataset[dataset$block_number != 0,] practice_counter = dataset %>% group_by(PROLIFIC_PID) %>% summarise(counter = unique(counter_practice)) if (!suppressPlot){ windows() barplot(practice_counter$counter, names= as.character(1:length(practice_counter$counter))) title("Number of times practice phase") } print("Number of times practice phase summary:") print(summary(practice_counter$counter)) # check number of trials if (length(dataset$PROLIFIC_PID)/length(unique(dataset$PROLIFIC_PID)) == 72*10){ n_trials_per_subj = 72*10 print("Number of trials check --> ok") } else { stop("ERROR: number of trials are not 720 (72*10) 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 ################## # check for keyboard response lagging (time from mask to accepted response, shoud be 0) summary(as.numeric(dataset$time_keyboard_response) - as.numeric(dataset$time_mask)) # mask duration should be 0 # check for duration of target and delay (expected 60 + 100 ms) summary(as.numeric(dataset$time_delay_target) - as.numeric(dataset$time_target)) # duration of target presentation should be 60 summary(as.numeric(dataset$time_mask) - as.numeric(dataset$time_delay_target)) # duration of delay should be 100 # Add to the reaction times both the time of target presentation and the time of the delay 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) dataset$rt = dataset$reaction_time # ####### trial outliers ######### # dataset$rt_outliers_logic = IQR_custom("rt", "subj_nr", dataset) # #dataset$rt_outliers_logic[which(dataset$SOA_ms > 1000)] = FALSE # # # % of trials to be discarded as outliers in total # (sum(!dataset$rt_outliers_logic ))/length(dataset$rt_outliers_logic )*100 # # # % outliers per subj # dataset %>% # group_by(subj_nr) %>% # dplyr::summarise( # outliers_per_subj = sum(!rt_outliers_logic, na.rm=TRUE)*100/number_of_trials) %>% # as.data.frame() if (!suppressPlot){ #### show mean RT per block per subj ##### dataset$block_number = as.factor(dataset$block_number) windows() p = dataset%>% ggplot(aes(x=block_number, y=rt ))+ geom_boxplot()+ facet_grid(.~subj_nr )+ ggtitle("Mean RT per block per subj")+ theme(plot.title = element_text(hjust = 0.5))+ xlab("Block number") + ylab("RT") print(p) #### show accuracy per subj(block) ######## windows() p = dataset%>% group_by(block_number,subj_nr)%>% summarise(accuracy = mean(correct), subj_nr)%>% ggplot(aes(x=subj_nr, y=accuracy))+ geom_boxplot()+ # facet_grid(.~subj_nr, labeller = labeller(.cols = label_both))+ ggtitle("Accuracy per participant\n(one value per participant)")+ theme(plot.title = element_text(hjust = 0.5))+ xlab("Participant number") + ylab("Accuracy") print(p) #### show accuracy per block (subj) ######## windows() p = dataset%>% group_by(block_number,subj_nr)%>% summarise(accuracy = mean(correct), block_number)%>% ggplot(aes(x=block_number, y=accuracy))+ geom_boxplot()+ # facet_grid(.~subj_nr, labeller = labeller(.cols = label_both))+ ggtitle("Accuracy per block\n(one value per block)")+ theme(plot.title = element_text(hjust = 0.5))+ xlab("Block number") + ylab("Accuracy") print(p) } # Keep this for later (we need it to calculate inv_eff) # discard incorrect answers # dataset = dataset[as.logical(dataset$correct),] # & dataset$rt_outliers_logic,] # false are outliers # 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' dataset$Phase = dataset$block_number == 1 | dataset$block_number == 2 dataset$Phase[dataset$Phase] = "Bonus" dataset$Phase[dataset$Phase != "Bonus"] = "Earn" dataset$Phase = as.factor(dataset$Phase) return(dataset) } 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") ) } # prepare data dataset = prepare_data("exp3_dataset.csv", suppressPlot=TRUE) number_of_trials = 72*10 # declare number of trials in total number_of_subjs = length(unique(dataset$PROLIFIC_PID)) # number of subjects in dataset # check the % of trials over 1700ms reaction times by subj to check it's marginal dataset %>% # per target pos and subj group_by(subj_nr, target_pos) %>% summarise(how_many = sum(rt>1700)*100/length(rt)) %>% print(n=100) dataset %>% # per subj group_by(subj_nr) %>% summarise(how_many = sum(rt>1700)*100/length(rt)) %>% print(n=100) dataset %>% # general and sd group_by(subj_nr) %>% summarise(how_many = sum(rt>1700)*100/length(rt)) %>% summarise(mean(how_many), sd(how_many)) sum(dataset$rt>1700)*100/length(dataset$rt) # % trials over 1700ms dataset = dataset[dataset$rt<1700,] windows() plot_pp_outliers() dataset = dataset[!dataset$subj_nr %in% c(47,3,44),] dataset$subj_nr = factor(dataset$subj_nr) length(levels(dataset$subj_nr)) # Accuracy acc = dataset %>% group_by(subj_nr)%>% summarise(accuracy = mean(correct)) summary(acc$accuracy) #### plots obj-based effect #### # object based effect by phase {# object based-effect by Phase windows() dataset[as.logical(dataset$correct) & dataset$target_pos != "valid",] %>% group_by(target_pos, subj_nr, Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj_nr,Phase) %>% mutate(subj_mean = mean(rt), target_pos, n) %>% group_by(target_pos,Phase) %>% mutate(group_mean = mean(rt), subj_nr, n) %>% group_by(Phase) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(target_pos,Phase)%>% 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(600,800))+ theme_minimal()+ facet_grid(. ~ Phase)+ scale_x_discrete(labels= c("Different obj", "Same obj"))+ 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") } # object based-effect by block_number { windows() dataset[as.logical(dataset$correct) & dataset$target_pos != "valid",] %>% group_by(target_pos, subj_nr, block_number) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj_nr,block_number) %>% mutate(subj_mean = mean(rt), target_pos, n) %>% group_by(target_pos,block_number) %>% mutate(group_mean = mean(rt), subj_nr, n) %>% group_by(block_number) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(target_pos,block_number)%>% 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(600,820))+ theme_minimal()+ facet_grid(. ~ block_number)+ scale_x_discrete(labels= c("Different obj", "Same obj"))+ 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") } ############### INVALID TRIALS ################# #### 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) #### bar plots ###### #### Accuracy #### sum(obj_eff_rwd_data$correct)/length(obj_eff_rwd_data$correct) {windows() p = obj_eff_rwd_data %>% group_by(rwd, trg,subj, Phase) %>% summarise(correct = mean(correct), n=n()) %>% group_by(subj, Phase) %>% mutate(subj_mean = mean(correct), rwd, trg, n) %>% group_by(rwd, trg, Phase) %>% mutate(group_mean = mean(correct), subj, n) %>% group_by() %>% mutate(isv = (correct - subj_mean + group_mean)) %>% group_by(rwd, trg, Phase)%>% summarise(correct = mean(correct), isv = sem_fun(isv)) %>% ggplot(aes(x=rwd, y=correct, fill=trg))+ geom_bar(stat="identity", width=0.5, position = "dodge")+ geom_errorbar(aes(ymin=correct-isv, ymax=correct+isv), width=.2,position=position_dodge(.5))+ coord_cartesian(ylim = c(0.75,1))+ scale_y_continuous(breaks=seq(0.5,1, 0.1))+ theme_minimal()+ facet_grid(. ~ Phase)+ 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 interacting with reward\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("ACC") p } #### plot obj-based eff and reward #### # Earn phase { windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 & block_number != 2) %>% # earn phase 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(group_mean = mean(rt), subj, n,rwd, trg) %>% # this is an alternative, I think a wrong one but have to check (look at the top of the script for the link) 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,680))+ 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") p } # Phase wise with more breaks {windows() p = 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(subj, Phase) %>% mutate(subj_mean = mean(rt), rwd, trg, n) %>% group_by(rwd, trg, Phase) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, trg, Phase)%>% 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,870))+ scale_y_continuous(breaks=seq(600, 960, 20))+ theme_minimal()+ facet_grid(. ~ Phase)+ 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 interacting eith reward\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 } ##### violin plot ###### ## new plot { plot_data = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% group_by(rwd, trg,subj, Phase) %>% 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,Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(rwd, trg,Phase) %>% 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))+ facet_grid(. ~ Phase)+ coord_cartesian(ylim = c(400,1400))+ scale_y_continuous(breaks=seq(400, 1400, 100)) p } #### bar plots again ###### # Bonus phase { windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number == 1 | block_number == 2) %>% # Bonus phase 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(group_mean = mean(rt), subj, n,rwd, trg) %>% # this is an alternative, I think a wrong one but have to check (look at the top of the script for the link) 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(650,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") p } # Phase wise {windows() p = 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(subj, Phase) %>% mutate(subj_mean = mean(rt), rwd, trg, n) %>% group_by(rwd, trg, Phase) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, trg, Phase)%>% 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,900))+ theme_minimal()+ facet_grid(. ~ Phase)+ 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 interacting eith reward\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 } # Reward and phase only {windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% group_by(rwd, subj, Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj, Phase) %>% mutate(subj_mean = mean(rt), rwd, n) %>% group_by(rwd, Phase) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by(Phase) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, Phase)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=Phase, 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(600,900))+ theme_minimal()+ scale_x_discrete(labels= c("bonus", "extinction"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Low reward", "congruent", "incongruent"))+ ggtitle(paste("Object-based effect interacting eith reward\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 } # Reward and phase only BONUS {windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number == 1 | block_number == 2) %>% group_by(rwd, subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), rwd, n) %>% group_by(rwd) %>% mutate(group_mean = mean(rt), subj, n) %>% 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(600,900))+ theme_minimal()+ # scale_x_discrete(labels= c("bonus", "extinction"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Low reward", "congruent", "incongruent"))+ ggtitle(paste("Reward and phase only BONUS\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+ theme(plot.title = element_text(hjust = 0.5)) # xlab("Reward & Target") + ylab("RT") p } # Reward and phase only Extinction {windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 | block_number != 2) %>% group_by(rwd, subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), rwd, n) %>% group_by(rwd) %>% mutate(group_mean = mean(rt), subj, n) %>% 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(640,700))+ theme_minimal()+ # scale_x_discrete(labels= c("bonus", "extinction"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Low reward", "congruent", "incongruent"))+ ggtitle(paste("# Reward and phase only Extinction\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+ theme(plot.title = element_text(hjust = 0.5)) # xlab("Reward & Target") + ylab("RT") p } # target only Extinction {windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 | block_number != 2) %>% group_by(trg, subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), trg, n) %>% group_by(trg) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(trg)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=trg, 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(640,700))+ theme_minimal()+ # scale_x_discrete(labels= c("bonus", "extinction"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Low reward", "congruent", "incongruent"))+ ggtitle(paste("# Reward and phase only Extinction\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 - 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(block_number) %>% 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_continuous(breaks = c(1:10), minor_breaks = c(1:10))+ scale_fill_discrete(name = "Cue & Target", labels = c("Different obj", "Same obj"))+ ggtitle(paste("Object-based effect during the 10 blocks reward\n(error: intra-subj variability)\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+ theme(plot.title = element_text(hjust = 0.5))+ facet_grid(. ~ rwd)+ xlab("Reward & Target") + ylab("RT") p } # block number - line plot, all rwd conditions in the same window { 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 } #### inspect interaction in Bonus Phase #### # excluding low reward first obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number == 1 | block_number == 2) %>% # only Bonus Phase (block 1 & 2) subset(side_high_reward != "to_be_assigned") %>% # only High reward group_by(rwd, trg,subj) %>% summarise(rt = mean(rt)) %>% group_by(subj, trg) %>% summarize(rt = rt[rwd == "incongruent"] - rt[rwd == "congruent"]) %>% group_by(subj) %>% summarize(rt = rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>% summarize(mean(rt)) # plot objectXReward interaction index per subj Bonus library(forcats) windows() obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number == 1 | block_number == 2) %>% # only Bonus Phase (block 1 & 2) subset(side_high_reward != "to_be_assigned") %>% # only High reward group_by(rwd, trg,subj) %>% summarise(rt = mean(rt)) %>% group_by(subj, trg) %>% summarize(rt = rt[rwd == "incongruent"] - rt[rwd == "congruent"]) %>% group_by(subj) %>% summarize(rt = rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>% ggplot(aes(x=fct_reorder(subj,rt), y=rt))+ 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(640,700))+ theme_minimal() #### inspect interaction in Extinction Phase #### # excluding low reward first # plot objectXReward interaction index per subj Extinction library(forcats) windows() obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 | block_number != 2) %>% # only Bonus Phase (block 1 & 2) subset(side_high_reward != "to_be_assigned") %>% # only High reward group_by(rwd, trg,subj) %>% summarise(rt = mean(rt)) %>% group_by(subj, trg) %>% summarize(rt = rt[rwd == "incongruent"] - rt[rwd == "congruent"]) %>% group_by(subj) %>% summarize(rt = rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>% ggplot(aes(x=fct_reorder(subj,rt), y=rt))+ 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(640,700))+ theme_minimal() # plot objectXReward interaction index by block library(forcats) windows() obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(side_high_reward != "to_be_assigned") %>% # only High reward group_by(rwd, trg,subj,block_number) %>% summarise(rt = mean(rt)) %>% group_by(subj, trg, block_number) %>% summarize(rt = rt[rwd == "incongruent"] - rt[rwd == "congruent"]) %>% group_by(subj,block_number) %>% summarize(rt = rt[trg == "diff_obj"] - rt[trg == "same_obj"]) %>% group_by(block_number) %>% summarize(inter_indx = mean(rt))%>% ggplot(aes(x=block_number, y=inter_indx))+ 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(640,700))+ theme_minimal() subset(block_number == 1 | block_number == 2) %>% # only Bonus Phase (block 1 & 2) ### plot variablility #### {windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 | block_number != 2) %>% group_by(trg, subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), trg, n) %>% group_by(trg) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(trg)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=trg, 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(640,700))+ theme_minimal()+ # scale_x_discrete(labels= c("bonus", "extinction"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Low reward", "congruent", "incongruent"))+ ggtitle(paste("# Reward and phase only Extinction\nNumber of participants",length(levels(obj_eff_rwd_data$subj))))+ theme(plot.title = element_text(hjust = 0.5)) # xlab("Reward & Target") + ylab("RT") p } # normalize everything # by block object effect and rwd (c'รจ anche sopra con il 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, trg,block_number) %>% 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(600,680))+ 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(. ~ block_number) p } # by block object effect { windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% group_by(block_number, trg,subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), block_number, trg, n) %>% group_by(block_number, trg) %>% mutate(group_mean = mean(rt), subj, n) %>% # group_by() %>% # mutate(group_mean = mean(rt), subj, n,rwd, trg) %>% # this is an alternative, I think a wrong one but have to check (look at the top of the script for the link) group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(block_number, trg)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=as.factor(block_number), 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,810))+ 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("Block number") + ylab("RT") p } # by block reward { windows() p = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% group_by(block_number, rwd,subj) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), block_number, rwd, n) %>% group_by(block_number, rwd) %>% mutate(group_mean = mean(rt), subj, n) %>% # group_by() %>% # mutate(group_mean = mean(rt), subj, n,rwd, trg) %>% # this is an alternative, I think a wrong one but have to check (look at the top of the script for the link) group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(block_number, rwd)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=as.factor(block_number), 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(600,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("Block number") + ylab("RT") p } ########### STATS INVALID TRIALS ############### library(BayesFactor) set.seed(1000) options(contrasts = c("contr.sum","contr.poly")) # 3X2X2 RM-ANOVA per_anova_earn = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$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, within = c(rwd, trg, Phase), effect.size = "pes") get_anova_table(res.aov) per_anova_earn = as.data.frame(per_anova_earn) bf = anovaBF(rt ~ rwd*trg*Phase + subj, data = per_anova_earn, whichRandom="subj", iterations = 100000) bf bf/bf plot(bf) # 3X2X2 RM-ANOVA per_anova_earn = obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct) & obj_eff_rwd_data$Phase != "Bonus",] %>% group_by(subj, rwd,trg, block_number)%>% summarise(rt = mean(rt)) #str(per_anova_earn) per_anova_earn = ungroup(per_anova_earn) str(per_anova_earn) per_anova_earn$block_number = as.factor(per_anova_earn$block_number) res.aov <- anova_test(data = per_anova_earn, dv = rt, wid = subj, within = c(rwd, trg, block_number), effect.size = "pes") get_anova_table(res.aov) #### earn phase ##### obj_eff_earn = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Earn",] obj_eff_earn$Phase = factor(obj_eff_earn$Phase) # relevel #### ACC #### #### 2X2 anova: rew-trg and cue-trg --- high reward ##### # without low reward baseline obj_eff_high_rwd = obj_eff_earn[obj_eff_earn$side_high_reward != "to_be_assigned",] obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd) # relevel with only high rwd trials library(lme4) library(lmerTest) ciaobelli = glmer(correct ~ rwd * trg + (1 | subj), data = obj_eff_high_rwd, family = binomial(link = "logit")) #The link-function describing the relationship between your linear model and the binomial distribution of your outcome r_int<- ranef(model)$subj$`(Intercept)` qqnorm(r_int) qqline(r_int) shapiro.test(r_int) options(contrasts = c("contr.sum","contr.poly")) summary(ciaobelli) # with ANOVA after boxcox transformation library(MASS) per_anova_earn = obj_eff_high_rwd %>% group_by(subj, rwd,trg)%>% summarise(correct = mean(correct)) #str(per_anova_earn) per_anova_earn = ungroup(per_anova_earn) str(per_anova_earn) b <- boxcox(lm(per_anova_earn$correct ~ 1)) # boxcox on ACC lambda <- b$x[which.max(b$y)] # Exact lambda lambda per_anova_earn$correct_ed <- (per_anova_earn$correct ^ lambda - 1) / lambda hist(per_anova_earn$correct_ed) hist(per_anova_earn$correct) 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) per_anova_earn %>% group_by(trg)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(rwd)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(rwd,trg)%>% summarise(correct = mean(correct)) %>% ggplot(aes(y=correct,x=rwd,fill=trg))+ geom_bar(stat="identity", width=0.5, position = "dodge") #### RT #### #### 2 paired samples t.test: cue-trg --- low reward #### ### only LOW REW obj_eff_low_rew = obj_eff_earn[obj_eff_earn$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 # inv efficiency inv_same = obj_eff_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "same_obj" & obj_eff_earn_meanSj$rwd == "low_reward"] inv_diff = obj_eff_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "diff_obj" & obj_eff_earn_meanSj$rwd == "low_reward"] t.test(inv_same, inv_diff, paired = T) # t test ttestBF(inv_same, inv_diff, paired = T) # BF ### ONLY congruent obj_eff_cong_rwd = obj_eff_earn[obj_eff_earn$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_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "same_obj" & obj_eff_earn_meanSj$rwd == "congruent"] cong_diff = obj_eff_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "diff_obj" & obj_eff_earn_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_earn[obj_eff_earn$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_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "same_obj" & obj_eff_earn_meanSj$rwd == "incongruent"] incong_diff = obj_eff_earn_meanSj$inv_eff[obj_eff_earn_meanSj$trg == "diff_obj" & obj_eff_earn_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_earn[obj_eff_earn$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 = 1000000) 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_earn_meanSj[obj_eff_earn_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_earn[obj_eff_earn$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_earn_meanSj[obj_eff_earn_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_earn[obj_eff_earn$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) ciao = lmer(rt ~ rwd*trg + (1|subj),data = per_anova_earn, REML = T) anova(ciao) summary(ciao) ciao_lm = lm(rt ~ rwd*trg + subj,data = per_anova_earn) summary(ciao_lm) 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) # look at residuals model_1 = lmer(rt ~ rwd * trg + (1|subj), data = obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] ) summary(model_1) windows() plot(model_1) anova(model_1) model_1 = lmer(rt ~ trg*rwd + (1|subj), data = per_anova_earn ) summary(model_1) library(performance) check_collinearity(model_1) # inv efficiency per_anova_earn = obj_eff_earn_meanSj[obj_eff_earn_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) #### 3X2 #### obj_eff_earn per_anova_earn = obj_eff_earn[as.logical(obj_eff_earn$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) #### Bonus phase #### obj_eff_bonus = obj_eff_rwd_data[obj_eff_rwd_data$Phase == "Bonus",] obj_eff_bonus$Phase = factor(obj_eff_bonus$Phase) # relevel #### ACC #### library(BayesFactor) set.seed(1000) #### 2 paired samples t.test: cue-trg --- low reward #### ### only LOW REW obj_eff_low_rew = obj_eff_bonus[obj_eff_bonus$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 %>% group_by(subj, rwd,trg)%>% summarise(correct = mean(correct)) #str(per_anova_earn) per_anova_earn = ungroup(per_anova_earn) str(per_anova_earn) # correct A_ttest_same = per_anova_earn$correct[per_anova_earn$trg == "same_obj"] A_ttest_diff = per_anova_earn$correct[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 #### 2X2 anova: rew-trg and cue-trg --- high reward ##### # without low reward baseline obj_eff_high_rwd = obj_eff_bonus[obj_eff_bonus$side_high_reward != "to_be_assigned",] obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd) # relevel with only high rwd trials library(lme4) library(lmerTest) ciaobelli = glmer(correct ~ rwd * trg + (1 | subj), data = obj_eff_high_rwd, family = binomial(link = "logit")) #The link-function describing the relationship between your linear model and the binomial distribution of your outcome r_int<- ranef(model)$subj$`(Intercept)` qqnorm(r_int) qqline(r_int) shapiro.test(r_int) options(contrasts = c("contr.sum","contr.poly")) summary(ciaobelli) # with ANOVA after boxcox transformation library(MASS) per_anova_earn = obj_eff_high_rwd %>% group_by(subj, rwd,trg)%>% summarise(correct = mean(correct)) #str(per_anova_earn) per_anova_earn = ungroup(per_anova_earn) str(per_anova_earn) b <- boxcox(lm(per_anova_earn$correct ~ 1)) # boxcox on ACC lambda <- b$x[which.max(b$y)] # Exact lambda lambda per_anova_earn$correct_ed <- (per_anova_earn$correct ^ lambda - 1) / lambda hist(per_anova_earn$correct_ed) hist(per_anova_earn$correct) 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) per_anova_earn %>% group_by(trg)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(rwd)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(rwd,trg)%>% summarise(correct = mean(correct)) %>% ggplot(aes(y=correct,x=rwd,fill=trg))+ geom_bar(stat="identity", width=0.5, position = "dodge") ########## RT ######## library(BayesFactor) set.seed(1000) #### 2 paired samples t.test: cue-trg --- low reward #### ### only LOW REW obj_eff_low_rew = obj_eff_bonus[obj_eff_bonus$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 # inv efficiency inv_same = obj_eff_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "same_obj" & obj_eff_bonus_meanSj$rwd == "low_reward"] inv_diff = obj_eff_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "diff_obj" & obj_eff_bonus_meanSj$rwd == "low_reward"] t.test(inv_same, inv_diff, paired = T) # t test ttestBF(inv_same, inv_diff, paired = T) # BF ### ONLY congruent obj_eff_cong_rwd = obj_eff_bonus[obj_eff_bonus$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_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "same_obj" & obj_eff_bonus_meanSj$rwd == "congruent"] cong_diff = obj_eff_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "diff_obj" & obj_eff_bonus_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_bonus[obj_eff_bonus$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_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "same_obj" & obj_eff_bonus_meanSj$rwd == "incongruent"] incong_diff = obj_eff_bonus_meanSj$inv_eff[obj_eff_bonus_meanSj$trg == "diff_obj" & obj_eff_bonus_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_bonus[obj_eff_bonus$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 = 100000) windows() plot(bf) against_best = bf/bf plot(against_best[,3]) title("2X2\ncongruent vs baseline\nanovaBF", adj = 0) bf # inv efficiency per_anova_earn = obj_eff_bonus_meanSj[obj_eff_bonus_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_bonus[obj_eff_bonus$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") 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_bonus_meanSj[obj_eff_bonus_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_bonus[obj_eff_bonus$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 = 500000) 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) # look at residuals model_1 = lmer(rt ~ rwd * trg + (1|subj), data = obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] ) summary(model_1) windows() plot(model_1) # inv efficiency per_anova_earn = obj_eff_bonus_meanSj[obj_eff_bonus_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) #### 3X2 #### obj_eff_bonus 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) bf = anovaBF(rt ~ rwd*trg + subj, data = per_anova_earn, whichRandom="subj") bf 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) # object effect during extinction obj_eff_rwd_data[as.logical(obj_eff_rwd_data$correct),] %>% subset(block_number != 1 | block_number != 2) %>% group_by(rwd, subj, PROLIFIC_PID) %>% summarise(rt = mean(rt)) %>% group_by(subj, PROLIFIC_PID) %>% summarize(rt = rt[rwd == "congruent"] - rt[rwd == "incongruent"]) #### ACC 2X2X2 anova: rew-trg and cue-trg and phase --- 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 library(lme4) library(lmerTest) ciaobelli = glmer(correct ~ rwd * trg * Phase + (1 | subj), data = obj_eff_high_rwd, family = binomial(link = "logit")) #The link-function describing the relationship between your linear model and the binomial distribution of your outcome r_int<- ranef(model)$subj$`(Intercept)` qqnorm(r_int) qqline(r_int) shapiro.test(r_int) options(contrasts = c("contr.sum","contr.poly")) summary(ciaobelli) # with ANOVA after boxcox transformation library(MASS) per_anova_earn = obj_eff_high_rwd %>% group_by(subj, rwd,trg,Phase)%>% summarise(correct = mean(correct)) #str(per_anova_earn) per_anova_earn = ungroup(per_anova_earn) str(per_anova_earn) b <- boxcox(lm(per_anova_earn$correct ~ 1)) # boxcox on ACC lambda <- b$x[which.max(b$y)] # Exact lambda lambda per_anova_earn$correct_ed <- (per_anova_earn$correct ^ lambda - 1) / lambda hist(per_anova_earn$correct_ed) hist(per_anova_earn$correct) res.aov <- anova_test(data = per_anova_earn, dv = correct_ed, wid = subj, within = c(rwd, trg,Phase), effect.size = "pes") get_anova_table(res.aov) per_anova_earn %>% group_by(rwd,trg,Phase)%>% summarise(correct = mean(correct)) %>% ggplot(aes(y=correct,x=rwd,fill=trg))+ geom_bar(stat="identity", width=0.5, position = "dodge")+ facet_grid(. ~ Phase) per_anova_earn %>% group_by(rwd)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(Phase)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(trg)%>% summarise(correct = mean(correct)) per_anova_earn %>% group_by(rwd,Phase)%>% summarise(correct = mean(correct)) %>% ungroup()%>% summarise( a = correct[rwd == 'congruent' & Phase == "Bonus"]- correct[rwd == 'incongruent' & Phase == "Bonus"], b = correct[rwd == 'congruent' & Phase == "Earn"]- correct[rwd == 'incongruent' & Phase == "Earn"] ) #### RT 2X2X2 anova: rew-trg and cue-trg and phase --- 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,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, within = c(rwd, trg,Phase), effect.size = "pes") get_anova_table(res.aov) library(BayesFactor) bf = anovaBF(rt ~ rwd*trg*Phase + subj, data = per_anova_earn, whichRandom="subj", iterations = 100000) windows() plot(bf[c(11,15)]) against_best = bf/bf plot(against_best[,3]) bf against_best[11,16] residuals(attr(res.aov, "args")$model) # look at residuals model_1 = lmer(rt ~ rwd * trg * Phase + (1|subj), data = obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] ) summary(model_1) windows() plot(model_1) #model_1 = lmer(rt ~ rwd * trg * Phase + (1|subj)+(1+subj|rwd), data = obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] ) #summary(model_1) ## first two blocks of extinction # 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) & obj_eff_high_rwd$block_number %in% c(1,2,3,4),] %>% 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, within = c(rwd, trg,Phase), effect.size = "pes") get_anova_table(res.aov) # plot by block object and reward seperatly and only on high rewarded trials # 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 # reward by block_number { windows() obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] %>% group_by(rwd, subj, block_number) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj,block_number) %>% mutate(subj_mean = mean(rt), rwd, n) %>% group_by(rwd,block_number) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by(block_number) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd,block_number)%>% 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(600,820))+ theme_minimal()+ facet_grid(. ~ block_number)+ # scale_x_discrete(labels= c("Different obj", "Same obj"))+ # 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") } # object by block_number { windows() obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] %>% group_by(trg, subj, block_number) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj,block_number) %>% mutate(subj_mean = mean(rt), trg, n) %>% group_by(trg,block_number) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by(block_number) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(trg,block_number)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=trg, 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,820))+ theme_minimal()+ facet_grid(. ~ block_number)+ # scale_x_discrete(labels= c("Different obj", "Same obj"))+ # 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") } # object by block_number + facet for reward { windows() obj_eff_high_rwd[as.logical(obj_eff_high_rwd$correct),] %>% group_by(trg, subj, block_number, rwd) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj,block_number) %>% mutate(subj_mean = mean(rt),rwd, trg, n) %>% group_by(trg,block_number,rwd) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by(block_number) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(trg,block_number,rwd)%>% summarise(rt = mean(rt), isv = sem_fun(isv)) %>% ggplot(aes(x=trg, 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,820))+ theme_minimal()+ facet_grid(rwd ~ block_number)+ # scale_x_discrete(labels= c("Different obj", "Same obj"))+ # 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") } #### Dealing with heteroscedaticity in: 2X2X2 anova: rew-trg and cue-trg and phase --- high reward #### #### subtraction #### # without low reward baseline & correct obj_eff_high_rwd = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward != "to_be_assigned" & as.logical(obj_eff_rwd_data$correct),] obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd) # relevel with only high rwd trials high = obj_eff_high_rwd%>% group_by(subj, rwd,trg,Phase)%>% summarise(rt = mean(rt)) high = ungroup(high) str(high) # only low reward baseline & correct obj_eff_low_rew = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward == "to_be_assigned" & as.logical(obj_eff_rwd_data$correct),] obj_eff_low_rew$rwd = factor(obj_eff_low_rew$rwd) # relevel low = obj_eff_low_rew%>% group_by(subj, rwd,trg,Phase)%>% summarise(rt = mean(rt)) low = ungroup(low) str(low) high$rt[high$trg == "diff_obj" & high$Phase == "Bonus"] = high$rt[high$trg == "diff_obj" & high$Phase == "Bonus"]-low$rt[low$trg == "diff_obj" & low$Phase == "Bonus"] high$rt[high$trg == "same_obj" & high$Phase == "Bonus"] = high$rt[high$trg == "same_obj" & high$Phase == "Bonus"]-low$rt[low$trg == "same_obj" & low$Phase == "Bonus"] high$rt[high$trg == "diff_obj" & high$Phase == "Earn"] = high$rt[high$trg == "diff_obj" & high$Phase == "Earn"]-low$rt[low$trg == "diff_obj" & low$Phase == "Earn"] high$rt[high$trg == "same_obj" & high$Phase == "Earn"] = high$rt[high$trg == "same_obj" & high$Phase == "Earn"]-low$rt[low$trg == "same_obj" & low$Phase == "Earn"] ## Statistical Model res.aov <- anova_test(data = high, dv = rt, wid = subj, within = c(rwd, trg,Phase), effect.size = "pes") get_anova_table(res.aov) # Phase wise {windows() p = high %>% group_by(rwd, trg,subj, Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj, Phase) %>% mutate(subj_mean = mean(rt), rwd, trg, n) %>% group_by(rwd, trg, Phase) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, trg, Phase)%>% 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,900))+ theme_minimal()+ facet_grid(. ~ Phase)+ 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 interacting eith reward\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 } #### division #### # without low reward baseline & correct obj_eff_high_rwd = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward != "to_be_assigned" & as.logical(obj_eff_rwd_data$correct),] obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd) # relevel with only high rwd trials high = obj_eff_high_rwd%>% group_by(subj, rwd,trg,Phase)%>% summarise(rt = mean(rt)) high = ungroup(high) str(high) # only low reward baseline & correct obj_eff_low_rew = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward == "to_be_assigned" & as.logical(obj_eff_rwd_data$correct),] obj_eff_low_rew$rwd = factor(obj_eff_low_rew$rwd) # relevel low = obj_eff_low_rew%>% group_by(subj, rwd,trg,Phase)%>% summarise(rt = mean(rt)) low = ungroup(low) str(low) high$rt[high$trg == "diff_obj" & high$Phase == "Bonus"] = high$rt[high$trg == "diff_obj" & high$Phase == "Bonus"]/((low$rt[low$trg == "diff_obj" & low$Phase == "Bonus"]+low$rt[low$trg == "same_obj" & low$Phase == "Bonus"])/2) high$rt[high$trg == "same_obj" & high$Phase == "Bonus"] = high$rt[high$trg == "same_obj" & high$Phase == "Bonus"]/((low$rt[low$trg == "diff_obj" & low$Phase == "Bonus"]+low$rt[low$trg == "same_obj" & low$Phase == "Bonus"])/2) high$rt[high$trg == "diff_obj" & high$Phase == "Earn"] = high$rt[high$trg == "diff_obj" & high$Phase == "Earn"]/((low$rt[low$trg == "diff_obj" & low$Phase == "Earn"]+low$rt[low$trg == "same_obj" & low$Phase == "Earn"])/2) high$rt[high$trg == "same_obj" & high$Phase == "Earn"] = high$rt[high$trg == "same_obj" & high$Phase == "Earn"]/((low$rt[low$trg == "diff_obj" & low$Phase == "Earn"]+low$rt[low$trg == "same_obj" & low$Phase == "Earn"])/2) ## Statistical Model res.aov <- anova_test(data = high, dv = rt, wid = subj, within = c(rwd, trg,Phase), effect.size = "pes") get_anova_table(res.aov) # look at residuals model_1 = lmer(rt ~ rwd * trg * Phase + (1|subj), data = high ) summary(model_1) windows() plot(model_1) # Phase wise {windows() p = high %>% group_by(rwd, trg,subj, Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj, Phase) %>% mutate(subj_mean = mean(rt), rwd, trg, n) %>% group_by(rwd, trg, Phase) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, trg, Phase)%>% 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,900))+ theme_minimal()+ facet_grid(. ~ Phase)+ 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 interacting eith reward\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 } #### Exp 2 & 3 comparison ###### load("exp2_dataset.R") exp2_b = exp2[!exp2$block_number %in% c(1,2),] exp2_b exp2_b_r = exp2_b[exp2_b$side_high_reward != "to_be_assigned" & as.logical(exp2_b$correct),] exp2_b_r$rwd = factor(exp2_b_r$rwd) # relevel with only high rwd trials obj_eff_high_rwd = obj_eff_rwd_data[obj_eff_rwd_data$side_high_reward != "to_be_assigned" & as.logical(obj_eff_rwd_data$correct) & obj_eff_rwd_data$Phase == "Earn",] obj_eff_high_rwd$rwd = factor(obj_eff_high_rwd$rwd) # relevel with only high rwd trials exp2_b_r_1 = exp2_b_r %>% group_by(rwd,trg,subj)%>% summarise(rt = mean(rt))%>% as.data.frame() exp2_b_r_1$exp = rep(2,length(exp2_b_r_1[,1])) obj_eff_high_rwd_1 = obj_eff_high_rwd %>% group_by(rwd,trg,subj)%>% summarise(rt = mean(rt))%>% as.data.frame() obj_eff_high_rwd_1$exp = rep(3,length(obj_eff_high_rwd_1[,1])) between_exp = rbind2(exp2_b_r_1,obj_eff_high_rwd_1) between_exp$exp = factor(between_exp$exp) str(between_exp) ## Statistical Model trhee.way <- aov(rt ~ rwd * trg * exp, data = between_exp) summary(trhee.way) # Phase wise {windows() p = between_exp %>% group_by(rwd, trg,subj, exp) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj, exp) %>% mutate(subj_mean = mean(rt), rwd, trg, n) %>% group_by(rwd, trg, exp) %>% mutate(group_mean = mean(rt), subj, n) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, trg, exp)%>% 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,900))+ theme_minimal()+ facet_grid(. ~ exp)+ 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 interacting eith reward\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 } ############### VALID TRIALS ################## #### plots obj-based effect #### # object based effect by phase {# object based-effect by Phase windows() dataset[as.logical(dataset$correct),] %>% group_by(target_pos, subj_nr, Phase) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj_nr,Phase) %>% mutate(subj_mean = mean(rt), target_pos, n) %>% group_by(target_pos,Phase) %>% mutate(group_mean = mean(rt), subj_nr, n) %>% group_by(Phase) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(target_pos,Phase)%>% 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(600,800))+ theme_minimal()+ facet_grid(. ~ Phase)+ scale_x_discrete(labels= c("Different obj", "Same obj", "Valid"))+ scale_fill_discrete(name = "Type trial", labels = c("Different obj", "Same obj", "Valid"))+ ggtitle("Object-based effect\n(error bar: intra-subj variability)")+ theme(plot.title = element_text(hjust = 0.5))+ xlab("Trial type") + ylab("RT") } # object based-effect by block_number { windows() dataset[as.logical(dataset$correct),] %>% group_by(target_pos, subj_nr, block_number) %>% summarise(rt = mean(rt), n=n()) %>% group_by(subj_nr,block_number) %>% mutate(subj_mean = mean(rt), target_pos, n) %>% group_by(target_pos,block_number) %>% mutate(group_mean = mean(rt), subj_nr, n) %>% group_by(block_number) %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(target_pos,block_number)%>% 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(600,820))+ theme_minimal()+ facet_grid(. ~ block_number)+ scale_x_discrete(labels= c("Different obj", "Same obj", "Valid"))+ scale_fill_discrete(name = "Type trial", labels = c("Different obj", "Same obj", "Valid"))+ ggtitle("Object-based effect\n(error bar: intra-subj variability)")+ theme(plot.title = element_text(hjust = 0.5))+ xlab("Trial type") + ylab("RT") } valid_rwd_data = dataset[dataset$target_pos == "valid",] # taking into account only valid trials valid_rwd_data$target_pos = as.factor(as.character(valid_rwd_data$target_pos)) #### added to trim invalid 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) # plot RT, rew-trg and cue-trg #### { windows() p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% group_by(rwd, subj, Phase) %>% summarise(rt = mean(rt)) %>% group_by(subj) %>% mutate(subj_mean = mean(rt), rwd, Phase) %>% group_by(rwd, Phase) %>% mutate(group_mean = mean(rt), subj) %>% group_by() %>% mutate(isv = (rt - subj_mean + group_mean)) %>% group_by(rwd, Phase)%>% 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(600,850))+ 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(600, 850, 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(. ~ Phase) p } ################ STATS VALID TRIALS #################### ### Phase x rwd #### per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% group_by(subj,rwd,Phase)%>% summarise(rt = mean(rt)) valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% group_by(subj,rwd,Phase)%>% summarise(rt = mean(rt)) %>% group_by(Phase,rwd) %>% summarise(rt = mean(rt)) valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% group_by(subj,rwd,Phase)%>% summarise(rt = mean(rt)) %>% group_by(subj,rwd) %>% summarise(rt = mean(rt))%>% group_by(rwd) %>% summarise(rt = mean(rt)) valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% group_by(subj,rwd,Phase)%>% summarise(rt = mean(rt)) %>% group_by(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, within = c(rwd,Phase), effect.size = "pes") get_anova_table(res.aov, correction = "GG") #### ACC #### { windows() p = valid_rwd_data %>% group_by(rwd, trg,subj, Phase) %>% summarise(acc = mean(correct)) %>% group_by(subj) %>% mutate(subj_mean = mean(acc), rwd, trg) %>% group_by(rwd, trg,Phase) %>% mutate(group_mean = mean(acc), subj) %>% group_by() %>% mutate(isv = (acc - subj_mean + group_mean)) %>% group_by(rwd, trg,Phase)%>% summarise(acc = mean(acc), isv = sem_fun(isv)) %>% ggplot(aes(x=rwd, y=acc, 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.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")+ facet_grid(. ~ Phase) p } per_anova_earn = valid_rwd_data %>% group_by(subj,rwd,Phase)%>% summarise(acc = mean(correct)) valid_rwd_data %>% group_by(subj,rwd,Phase)%>% summarise(acc = mean(correct))%>% group_by(Phase,rwd)%>% summarise(acc = mean(acc)) valid_rwd_data %>% group_by(subj,rwd,Phase)%>% summarise(acc = mean(correct))%>% group_by(Phase)%>% summarise(acc = mean(acc)) valid_rwd_data %>% group_by(subj,rwd,Phase)%>% summarise(acc = mean(correct))%>% group_by(rwd)%>% summarise(acc = mean(acc)) #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, Phase), effect.size = "pes") get_anova_table(res.aov, correction = "GG") #### Phase 1: rwd #### { windows() p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Bonus") %>% 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(600,850))+ 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(600, 850, 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 } per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Bonus")%>% group_by(subj,rwd)%>% summarise(rt = mean(rt)) valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Bonus")%>% group_by(subj,rwd)%>% summarise(rt = mean(rt))%>% group_by(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, correction = "GG") pwc <- per_anova_earn %>% pairwise_t_test( rt ~ rwd, paired = TRUE, p.adjust.method = "bonferroni", detailed = T ) data.frame(pwc[1,"p"]) 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) #### ACC #### { windows() p = valid_rwd_data %>% subset(Phase == "Bonus")%>% group_by(rwd, trg,subj) %>% summarise(acc = mean(correct)) %>% group_by(subj) %>% mutate(subj_mean = mean(acc), rwd, trg) %>% group_by(rwd, trg) %>% mutate(group_mean = mean(acc), subj) %>% group_by() %>% mutate(isv = (acc - subj_mean + group_mean)) %>% group_by(rwd, trg)%>% summarise(acc = mean(acc), isv = sem_fun(isv)) %>% ggplot(aes(x=rwd, y=acc, 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.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 } per_anova_earn = valid_rwd_data %>% subset(Phase == "Bonus")%>% group_by(subj,rwd)%>% summarise(acc = mean(correct)) valid_rwd_data %>% subset(Phase == "Bonus")%>% group_by(subj,rwd)%>% summarise(acc = mean(correct))%>% group_by(rwd)%>% summarise(acc = mean(acc)) #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) ### Phase 2: rwd ##### { windows() p = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Earn") %>% 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,650))+ 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, 650, 10))+ 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 } per_anova_earn = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Earn")%>% group_by(subj,rwd)%>% summarise(rt = mean(rt)) valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Earn")%>% group_by(subj,rwd)%>% summarise(rt = mean(rt))%>% group_by(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, correction = "GG") #### ACC #### { windows() p = valid_rwd_data %>% subset(Phase == "Earn")%>% group_by(rwd, trg,subj) %>% summarise(acc = mean(correct)) %>% group_by(subj) %>% mutate(subj_mean = mean(acc), rwd, trg) %>% group_by(rwd, trg) %>% mutate(group_mean = mean(acc), subj) %>% group_by() %>% mutate(isv = (acc - subj_mean + group_mean)) %>% group_by(rwd, trg)%>% summarise(acc = mean(acc), isv = sem_fun(isv)) %>% ggplot(aes(x=rwd, y=acc, 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.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 } per_anova_earn = valid_rwd_data %>% subset(Phase == "Earn")%>% group_by(subj,rwd)%>% summarise(acc = mean(correct)) valid_rwd_data %>% subset(Phase == "Earn")%>% group_by(subj,rwd)%>% summarise(acc = mean(correct))%>% group_by(rwd)%>% summarise(acc = mean(acc)) #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") ### inverse efficiency #### #### inv eff #### { windows() per_anova_earn %>% group_by(subj) %>% mutate(subj_mean = mean(inv_eff), rwd, n) %>% group_by(rwd) %>% mutate(group_mean = mean(inv_eff), subj, n) %>% group_by() %>% mutate(isv = (inv_eff - subj_mean + group_mean)) %>% group_by(rwd)%>% summarise(inv_eff = mean(inv_eff), isv = sem_fun(isv)) %>% ggplot(aes(y=inv_eff,x=rwd,fill=rwd))+ 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,700))+ scale_y_continuous(breaks=seq(600,700, 100))+ scale_x_discrete(labels= c("Baseline","Low reward", "High reward"))+ scale_fill_discrete(name = "Cue & Target", labels = c("Baseline","Low reward", "High reward"))+ ggtitle("Inverse efficiency score")+ geom_text(aes(label=round(inv_eff)), position=position_dodge(width=0.9), vjust=-0.25) } per_anova_earn = valid_rwd_data %>% subset(Phase == "Earn")%>% group_by(subj, rwd)%>% summarise(acc = mean(correct), n = n()) supporto = valid_rwd_data[as.logical(valid_rwd_data$correct),] %>% subset(Phase == "Earn")%>% group_by(subj, rwd)%>% 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) res.aov <- anova_test(data = per_anova_earn, dv = inv_eff, wid = subj, within = c(rwd), effect.size = "pes") get_anova_table(res.aov) per_anova_earn %>% group_by(rwd) %>% summarise(inv_eff = mean(inv_eff)) pwc <- per_anova_earn %>% pairwise_t_test( inv_eff ~ rwd, paired = TRUE, p.adjust.method = "bonferroni", detailed = T ) data.frame(pwc) library(effsize) cohen.d(per_anova_earn$inv_eff[per_anova_earn$rwd == "low_reward"],per_anova_earn$inv_eff[per_anova_earn$rwd == "congruent"], paired = T) cohen.d(per_anova_earn$inv_eff[per_anova_earn$rwd == "incongruent"],per_anova_earn$inv_eff[per_anova_earn$rwd == "congruent"], paired = T) cohen.d(per_anova_earn$inv_eff[per_anova_earn$rwd == "incongruent"],per_anova_earn$inv_eff[per_anova_earn$rwd == "low_reward"], paired = T)