# packages & workspace ------------------------------------------------------- library(tidyverse) library(lme4) library(emmeans) # clear workspace rm(list=ls()) # set working directory to wherever the r script is located setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) getwd() # load workspace load("R workspace ordinality pigs.RData") #load data from data set if not from workspace: op <- read.csv(file="Results ordinal pigs.csv", header=T, dec=",", sep=";") #save.image("R workspace ordinality pigs.RData") #op <- read.csv(file="C:/INSERT DIRECTORY/Results ordinal pigs.csv", header=T, dec=",", sep=";") summary(op) op <- op[1:960,] tail(op) op$Subject <- as.factor(op$Subject) op$Condition <- as.factor(op$Condition) op$Experimenter <- as.factor(op$Experimenter) op$Training_group <- as.factor(op$Training_group) op$Choice_number <- as.factor(op$Choice_number) op$Condition <- relevel(op$Condition, ref="Spatial") op$Condition <- relevel(op$Condition, ref="Sagittal") #one dataset per condition/test op_FP <- subset(op, op$Condition=="FP") op_sag <- subset(op, op$Condition=="Sagittal") op_spat <- subset(op, op$Condition=="Spatial") op_sag$Choice_number <- droplevels(op_sag$Choice_number) op_spat$Choice_number <- droplevels(op_spat$Choice_number) op_FP$Choice_number <- droplevels(op_FP$Choice_number) # paper plots ------------------------------------------------------------------ ## sagittal ---------------------- # to add raw data points to plot: # raw data converted to proportions per subject and bowl position op_sag_prop <- op_sag %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_sag_prop_exp <- expand.grid(Subject=unique(op_sag$Subject), Choice_number=c(1:6)) op_sag_prop_exp <- merge(op_sag_prop_exp, op_sag_prop, all.x = T) op_sag_prop_exp$n[is.na(op_sag_prop_exp$n)] <- 0 op_sag_prop_exp$Choice_prop[is.na(op_sag_prop_exp$Choice_prop)] <- 0 op_sag_prop_exp$Training_group <- ifelse(op_sag_prop_exp$Subject %in% unique(op_sag_prop$Subject[op_sag_prop$Training_group==5]), 5, 2) op_sag_prop_exp$Training_group <- as.factor(op_sag_prop_exp$Training_group) op_sag_prop_exp$Choice_number <- as.factor(op_sag_prop_exp$Choice_number) head(op_sag_prop_exp,10) str(op_sag_prop_exp) # df to colorize bowl positions df <-data.frame(x = factor(c(1,2,3,4,5,6)), colour = factor(c(1,2,3,4,5,1))) # plot 1a ggplot(data=op_sag_prop_exp, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + labs(title = "a) Sagittal", y = "Proportion of Bowl Choice", x = "Bowl Position", color="Group:") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F, show.legend = F) + scale_fill_manual(values=c("white", "mediumorchid4","white", "white", "deepskyblue3")) + geom_hline(yintercept=1/6, linetype="dashed") + #ylim(0,1) + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) # raw data converted to proportions per subject and bowl position op_sag_prop_first3 <- op_sag %>% filter(Trial<4) %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_sag_prop_exp_first3 <- expand.grid(Subject=unique(op_sag$Subject), Choice_number=c(1:6)) op_sag_prop_exp_first3 <- merge(op_sag_prop_exp_first3, op_sag_prop_first3, all.x = T) op_sag_prop_exp_first3$n[is.na(op_sag_prop_exp_first3$n)] <- 0 op_sag_prop_exp_first3$Choice_prop[is.na(op_sag_prop_exp_first3$Choice_prop)] <- 0 op_sag_prop_exp_first3$Training_group <- ifelse(op_sag_prop_exp_first3$Subject %in% unique(op_sag_prop_first3$Subject[op_sag_prop_first3$Training_group==5]), 5, 2) op_sag_prop_exp_first3$Training_group <- as.factor(op_sag_prop_exp_first3$Training_group) op_sag_prop_exp_first3$Choice_number <- as.factor(op_sag_prop_exp_first3$Choice_number) head(op_sag_prop_exp_first3,10) # df to colorize bowl positions df2 <-data.frame(x = factor(c(1,2,3,4,5,6)), colour = factor(c(1,2,3,4,5,1))) df2$colour <- factor(df2$colour, levels = c(1,2,5,3,4)) # plot 1b # op_firsts_sag_mean ggplot(data=op_sag_prop_exp_first3, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + labs(title = "b) Sagittal - First 3 Trials", y = "Proportion of Bowl Choice", x = "Bowl Position", fill="") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df2, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F) + scale_fill_manual(values=c("white", "mediumorchid4", "deepskyblue3","white", "white"), labels=c("", "", "Target Bowls", "", "")) + geom_hline(yintercept=1/6, linetype="dashed") + #ylim(0,1) + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "left", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) ## sagittal errors ----------------------- # plot of errors errors <- data.frame(x = c(1, 2, 3, 4, 5, 6), colour = factor(c(1,5,3,2,5,4))) op_sag %>% filter(Choice_correct==0) %>% group_by(Subject, Training_group) %>% count(Choice_number, .drop=F) %>% mutate(prop=prop.table(n)) %>% filter(!is.na(prop)) %>% mutate(n2=ifelse(as.character(Training_group)==as.character(Choice_number), NA, n)) %>% ggplot(aes(x = Choice_number, y = n2, color = Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) + geom_tile(data=errors, aes(x = x, y = -.6, fill = colour, height = .3), inherit.aes = F) + scale_fill_manual(values=c("mediumorchid4", "deepskyblue3","plum1","lightskyblue1", "white"), labels = c("", "before", "", "after", "")) + scale_y_continuous(expand = c(0,0), limits = c(-.9,8.5))+ labs(x="Bowl Position", y="Nr. of Errors", color="Group:", fill="Error:", title = "Sagittal Test") op_sag %>% filter(Choice_correct==0) %>% group_by(Training_group) %>% count(Choice_number, .drop=F) %>% mutate(prop=prop.table(n)) %>% filter(!is.na(prop)) ## spatial -------------------------- # to add raw data points to plot: # raw data converted to proportions per subject and bowl position op_spat_prop <- op_spat %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_spat_prop_exp <- expand.grid(Subject=unique(op_spat$Subject), Choice_number=c(1:6)) op_spat_prop_exp <- merge(op_spat_prop_exp, op_spat_prop, all.x = T) op_spat_prop_exp$n[is.na(op_spat_prop_exp$n)] <- 0 op_spat_prop_exp$Choice_prop[is.na(op_spat_prop_exp$Choice_prop)] <- 0 op_spat_prop_exp$Training_group <- ifelse(op_spat_prop_exp$Subject %in% unique(op_spat_prop$Subject[op_spat_prop$Training_group==5]), 5, 2) op_spat_prop_exp$Training_group <- as.factor(op_spat_prop_exp$Training_group) op_spat_prop_exp$Choice_number <- as.factor(op_spat_prop_exp$Choice_number) head(op_spat_prop_exp,10) str(op_spat_prop_exp) # df to colorize bowl positions df <-data.frame(x = factor(c(1,2,3,4,5,6)), colour = factor(c(1,2,3,4,5,1))) # plot 2a ggplot(data=op_spat_prop_exp, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F, show.legend = F) + scale_fill_manual(values=c("white", "mediumorchid4","plum", "lightskyblue", "deepskyblue3")) + labs(title = "a) Ordinal-Spatial Incongruence", y = "Proportion of Bowl Choice", x = "Bowl Position", color="Group:") + geom_hline(yintercept=1/6, linetype="dashed") + ylim(0,1) + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) # raw data converted to proportions per subject and bowl position op_spat_prop_first3 <- op_spat %>% filter(Trial<4) %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_spat_prop_exp_first3 <- expand.grid(Subject=unique(op_spat$Subject), Choice_number=c(1:6)) op_spat_prop_exp_first3 <- merge(op_spat_prop_exp_first3, op_spat_prop_first3, all.x = T) op_spat_prop_exp_first3$n[is.na(op_spat_prop_exp_first3$n)] <- 0 op_spat_prop_exp_first3$Choice_prop[is.na(op_spat_prop_exp_first3$Choice_prop)] <- 0 op_spat_prop_exp_first3$Training_group <- ifelse(op_spat_prop_exp_first3$Subject %in% unique(op_spat_prop_first3$Subject[op_spat_prop_first3$Training_group==5]), 5, 2) op_spat_prop_exp_first3$Training_group <- as.factor(op_spat_prop_exp_first3$Training_group) op_spat_prop_exp_first3$Choice_number <- as.factor(op_spat_prop_exp_first3$Choice_number) head(op_spat_prop_exp_first3,10) # plot 2b df2 <-data.frame(x = factor(c(1,2,3,4,5,6)), colour = factor(c(1,2,3,4,5,1))) df2$colour <- factor(df2$colour, levels = c(1,2,5,3,4)) # op_firsts_spat_mean ggplot(data=op_spat_prop_exp_first3, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + labs(title = "b) Ordinal-Spatial Inc. - First 3 Trials", y = "Proportion of Bowl Choice", x = "Bowl Position", fill="") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df2, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F) + scale_fill_manual(values=c("white", "mediumorchid4", "deepskyblue3","plum", "lightskyblue"), labels=c("", "", "Target Bowls ", "", "Spatial Distractor Bowls")) + geom_hline(yintercept=1/6, linetype="dashed") + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) ## fronto-parallel ---------------------- # to add raw data points to plot: # raw data converted to proportions per subject and bowl position op_FP_prop <- op_FP %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_FP_prop_exp <- expand.grid(Subject=unique(op_FP$Subject), Choice_number=levels(op_FP$Choice_number)[2:7]) op_FP_prop_exp <- merge(op_FP_prop_exp, op_FP_prop, all.x = T) op_FP_prop_exp$n[is.na(op_FP_prop_exp$n)] <- 0 op_FP_prop_exp$Choice_prop[is.na(op_FP_prop_exp$Choice_prop)] <- 0 op_FP_prop_exp$Training_group <- ifelse(op_FP_prop_exp$Subject %in% unique(op_FP_prop$Subject[op_FP_prop$Training_group==5]), 5, 2) op_FP_prop_exp$Training_group <- as.factor(op_FP_prop_exp$Training_group) op_FP_prop_exp$Choice_number <- factor(op_FP_prop_exp$Choice_number, labels = c("1L|6R", "2L|5R", "3L|4R", "4L|3R", "5L|2R", "6L|1R"), levels = c("1L", "2L", "3L", "3R", "2R", "1R")) head(op_FP_prop_exp,10) str(op_FP_prop_exp) # df to colorize bowl positions df3 <-data.frame(x = c(1, 2,2.5, 3, 4, 5,5.5, 6), colour = factor(c(1,2,3,1,1,4,5,1))) # plot 2a ggplot(data=op_FP_prop_exp, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df3, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F, show.legend = F) + scale_fill_manual(values=c("white", "mediumorchid4", "lightskyblue1", "deepskyblue3","plum1")) + labs(title = "a) Fronto-Parallel", y = "Proportion of Bowl Choice", x = "Bowl Position", color="Group:") + geom_hline(yintercept=1/3, linetype="dashed") + ylim(0,1) + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) # raw data converted to proportions per subject and bowl position op_FP_prop_first3 <- op_FP %>% filter(Trial<4) %>% group_by(Subject, Training_group) %>% count(Choice_number) %>% mutate(Choice_prop=prop.table(n)) # add 0% for non-chosen positions: op_FP_prop_exp_first3 <- expand.grid(Subject=unique(op_FP$Subject), Choice_number=levels(op_FP$Choice_number)[2:7]) op_FP_prop_exp_first3 <- merge(op_FP_prop_exp_first3, op_FP_prop_first3, all.x = T) op_FP_prop_exp_first3$n[is.na(op_FP_prop_exp_first3$n)] <- 0 op_FP_prop_exp_first3$Choice_prop[is.na(op_FP_prop_exp_first3$Choice_prop)] <- 0 op_FP_prop_exp_first3$Training_group <- ifelse(op_FP_prop_exp_first3$Subject %in% unique(op_FP_prop_first3$Subject[op_FP_prop_first3$Training_group==5]), 5, 2) op_FP_prop_exp_first3$Training_group <- as.factor(op_FP_prop_exp_first3$Training_group) op_FP_prop_exp_first3$Choice_number <- factor(op_FP_prop_exp_first3$Choice_number, levels = c("1L", "2L", "3L", "3R", "2R", "1R"), labels = c("1L|6R", "2L|5R", "3L|4R", "4L|3R", "5L|2R", "6L|1R")) head(op_FP_prop_exp_first3,10) df4 <-data.frame(x = c(1, 2,2.5, 3, 4, 5,5.5, 6), colour = factor(c(1,2,3,1,1,4,5,1))) df4$colour <- factor(df4$colour, levels = c(2,4,3,5,1)) # plot 2b # op_firsts_FP_mean ggplot(data=op_FP_prop_exp_first3, aes(x=Choice_number, y=Choice_prop, col=Training_group)) + geom_point(position = position_jitterdodge(dodge.width = .3), alpha=.2, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "point",fun=mean, size=3, shape=19, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "line",fun=mean, size=1, show.legend = F) + stat_summary(aes(group=Training_group),position=position_dodge(w=.3),geom = "errorbar",fun.data = mean_se, size=1, width=.3, show.legend = F) + labs(title = "b) Fronto-Parallel - First 3 Trials", y = "Proportion of Bowl Choice", x = "Bowl Position", fill="") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + geom_tile(data=df4, aes(x = x, y = -.06, fill = colour, height = .03), inherit.aes = F) + scale_fill_manual(values=c("mediumorchid4", "deepskyblue3","lightskyblue1", "plum1", "white"), labels=c("", "Left Search Bias ", "", "Right Search Bias", "")) + geom_hline(yintercept=1/3, linetype="dashed") + scale_y_continuous(expand = c(0,0), limits = c(-.09,1.05)) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "left", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) # Analysis 1 - training -------------------------------------------------------- str(op) op_train <- op %>% select(Subject, Training_group, Training_sessions, Training_trials) %>% unique.data.frame() op_train$sex <- c("m","m","f","f","m","m","f","f","m","f","m","m","f","f","m","m") # descriptives ................................................................. mean(op_train$Training_sessions) sd(op_train$Training_sessions) ## pigs needed on average 6.7 +/- 3.2 sessions mean(op_train$Training_trials/op_train$Training_sessions) sd(op_train$Training_trials/op_train$Training_sessions) ## pigs needed on average 19.1 +/- 2.9 trials per session mean(op_train$Training_trials) sd(op_train$Training_trials) ## pigs needed on average 129.4 +/- 68 trials # test if groups differ in their average number of trials ...................... # trial data is count data -> poisson full_train <- glm(Training_trials ~ Training_group, data = op_train, family = poisson) summary(full_train) # Overdispersion test AER::dispersiontest(full_train) # data: full_train # z = 4.5418, p-value = 2.789e-06 # alternative hypothesis: true dispersion is greater than 1 # sample estimates: # dispersion # 31.35293 ## data is over-dispersed # test goodness of fit 1 - pchisq(summary(full_train)$deviance, summary(full_train)$df.residual) # 0 ## The GOF test indicates that the Poisson model does not fit the data (p < 0.05) ## Use negbin model library(MASS) full_train.nb <- glm.nb(Training_trials ~ Training_group, data = op_train) # test goodness of fit 1 - pchisq(summary(full_train.nb)$deviance, summary(full_train.nb)$df.residual) # 0.2195781 ## The GOF test indicates that the Neg.Bin. model fits the data (p > 0.05) summary(full_train.nb) # Coefficients: # Estimate Std. Error z value Pr(>|z|) # (Intercept) 4.7658 0.1890 25.216 <2e-16 *** # Training_group5 0.1670 0.2517 0.664 0.507 ## no effect of training group # predict mean counts nd <- data.frame(Training_group=unique(op_train$Training_group)) full_train_res <- cbind(nd, Mean = unique(predict(full_train.nb, newdata = nd, type="response")), SE = unique(predict(full_train.nb, newdata = nd, type="response", se.fit = TRUE)$se.fit)) # Training_group Mean SE # 1 2 117.4286 22.19371 # 2 5 138.7778 23.07082 # model plot ggplot(full_train_res, aes(x = Training_group, y = Mean, colour = Training_group)) + # geom_boxplot(data = op_train, mapping = aes(x = Training_group, y = Training_trials, colour = Training_group), alpha = .3, show.legend = F) + geom_point(data = op_train, mapping = aes(x = Training_group, y = Training_trials, colour = Training_group), alpha = .3, position = position_jitter(width = .2), show.legend = F) + geom_point(size=4, show.legend = F) + geom_errorbar(aes(x = Training_group, ymin = Mean+SE, ymax = Mean-SE), show.legend = F, width = .25) + labs(x="Group", y="Nr. of Trials", title = "Training") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + theme_classic() + theme(axis.title = element_text(size=16), axis.text=element_text(size=12), title=element_text(size=16), legend.position="bottom", legend.justification = "right", legend.direction = "horizontal", axis.ticks = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), panel.grid.major.x = element_blank(), panel.grid.minor = element_blank()) # Analysis 2 - tests ----------------------------------------------------------- ## sagittal ----------------------------------------------------------- # data mdata_sag <- op_sag %>% dplyr::select(Subject, Training_group, Trial, Rewarded, Choice_correct, Choice_number) %>% mutate(Rewarded=as.factor(Rewarded)) %>% mutate(Rewarded.lag=lag(Rewarded)) %>% droplevels() mdata_sag$Rewarded.lag[mdata_sag$Trial==1] <- "Y" # scale covariate mdata_sag$Trial.z <- scale(mdata_sag$Trial) # dummy-code and center factor variable mdata_sag$Rewarded.lag.Y <- (as.numeric(mdata_sag$Rewarded.lag)-1)-mean(as.numeric(mdata_sag$Rewarded.lag)-1) mdata_sag$Rewarded.Y <- (as.numeric(mdata_sag$Rewarded)-1)-mean(as.numeric(mdata_sag$Rewarded)-1) str(mdata_sag) # 320 obs # full model full_sag0 <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z |Subject), data = mdata_sag, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_sag0)$varcor # remove correlations between random intercept and random slopes full_sag <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z || Subject), data = mdata_sag, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Rewarded + Trial.z, data = mdata_sag) car::vif(xx) max(car::vif(xx)) # 1.005038 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_sag$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_sag[mdata_sag$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_sag@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_sag)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_sag <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_sag,3) # orig min max # (Intercept) 1.666 1.580 1.814 # Training_group5 -1.049 -1.338 -0.697 # RewardedY -0.187 -0.299 -0.042 # Trial.z -0.313 -0.374 -0.254 ggplot(mod.stab.res_full_sag, aes(x=orig, y = row.names(mod.stab.res_full_sag)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## model is stable # inference: .......................................................... null_sag <- glmer(Choice_correct ~ #Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z || Subject), data = mdata_sag, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_sag, null_sag, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_sag 6 353.22 375.83 -170.61 341.22 # full_sag 7 352.77 379.15 -169.39 338.77 2.4472 1 0.1177 ## no significant effect of group on correct choice drop1(full_sag, test="Chisq") # npar AIC LRT Pr(Chi) # 352.77 # Training_group 1 353.22 2.4472 0.11774 # Rewarded 1 351.20 0.4323 0.51088 # Trial.z 1 355.41 4.6383 0.03127 * round(summary(full_sag)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) 1.666 0.494 # Training_group5 -1.049 0.634 # RewardedY -0.187 0.280 # Trial.z -0.313 0.141 # estimate and CI in response space per group emmeans(full_sag, pairwise ~ "Training_group", type="response") # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.828 0.0683 Inf 0.653 0.925 # 5 0.628 0.0975 Inf 0.427 0.793 ## in both groups the lower CI is above the chance level (0.16667) test(emmeans(full_sag, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.828 0.0683 Inf 0.167 6.628 <.0001 # 5 0.628 0.0975 Inf 0.167 5.109 <.0001 ## both groups perform significantly above chance in choosing the target bowl ### model plot ................................................................. mdata_sag_plot <- as.data.frame(effects::effect("Training_group", full_sag)) mdata_sag_plot %>% ggplot(aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "Sagittal Test", x="Group", y="Probability of choosing the target bowl") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ranef.data_sag <- as.data.frame(ggeffects::ggpredict(full_sag, terms = c("Trial.z[all]", "Training_group", "Subject"), type="random")) names(ranef.data_sag) <- c("Trial.z", "predicted", "std.error", "conf.low", "conf.high", "Training_group", "Subject") ranef.data_sag <- ranef.data_sag %>% filter(Training_group=="2" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="2"]) | Training_group=="5" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="5"])) %>% mutate(Subject=factor(Subject, levels = c(sort(unique(op_spat$Subject[op_spat$Training_group=="2"])), sort(unique(op_spat$Subject[op_spat$Training_group=="5"])))), unz.trial = round((Trial.z * attr(mdata_sag$Trial.z, 'scaled:scale') + attr(mdata_sag$Trial.z, 'scaled:center')),0)) ranef.data_sag %>% ggplot(aes(x = unz.trial, y = predicted))+ geom_hline(yintercept = 1/6, lty=2) + # chance level geom_line(aes(colour = Training_group), size=1) + # average probability geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Training_group), alpha=.2, show.legend = F) + geom_point(data=mdata_sag, aes(x = Trial, y = Choice_correct), inherit.aes = F, col="grey50") + facet_wrap(~Subject) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + scale_fill_manual(values=c("mediumorchid4","deepskyblue3")) + labs(title="Sagittal Test", x="Trials", y="Probability of choosing the target bowl", col="Group") + theme_bw(base_size = 15) ## sagittal - first ----------------------------------------------- # data mdata_sag_f <- op_sag %>% dplyr::select(Subject, Training_group, Trial, Choice_correct, Choice_number) %>% filter(Trial<4) %>% droplevels() # scale covariate mdata_sag_f$Trial.z <- scale(mdata_sag_f$Trial) str(mdata_sag_f) # 48 obs # full model full_sag_f0 <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_sag_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_sag_f0)$varcor # remove correlations between random intercept and random slopes full_sag_f <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z ||Subject), data = mdata_sag_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Trial.z, data = mdata_sag_f) car::vif(xx) max(car::vif(xx)) # 1 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_sag_f$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_sag_f[mdata_sag_f$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_sag_f@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(red.data, full.mod.stab, res, subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_sag_f)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_sag_f <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_sag_f,3) # orig min max # (Intercept) 2.586 2.275 3.310 # Training_group5 -1.327 -1.995 -0.646 # Trial.z -0.202 -0.499 0.000 ggplot(mod.stab.res_full_sag_f, aes(x=orig, y = row.names(mod.stab.res_full_sag_f)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## not super stable, but ok # inference: .......................................................... null_sag_f <- glmer(Choice_correct ~ #Training_group + Trial.z + (1 + Trial.z ||Subject), data = mdata_sag_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_sag_f, null_sag_f, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_sag_f 4 52.994 60.479 -22.497 44.994 # full_sag_f 5 53.259 62.615 -21.630 43.259 1.7346 1 0.1878 ## no significant effect of group on correct choice drop1(full_sag_f, test="Chisq") # npar AIC LRT Pr(Chi) # 53.259 # Training_group 1 52.994 1.73458 0.1878 # Trial.z 1 51.503 0.24351 0.6217 round(summary(full_sag_f)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) 2.586 0.996 # Training_group5 -1.327 1.052 # Trial.z -0.202 0.412 emmeans(full_sag_f, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.930 0.0648 Inf 0.654 0.989 # 5 0.779 0.1114 Inf 0.498 0.926 ## in both groups the lower CI is above the chance level (0.16667) test(emmeans(full_sag_f, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.930 0.0648 Inf 0.167 4.214 <.0001 # 5 0.779 0.1110 Inf 0.167 4.434 <.0001 ## both groups perform significantly above chance in choosing the target bowl # model plot ................................................................... mdata_sag_f_plot <- as.data.frame(effects::effect("Training_group", full_sag_f)) ggplot(mdata_sag_f_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "Sagittal Test - First 3 Trials", x="Group", y="Probability of choosing the target bowl") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ## spatial 2/5 -------------------------------------------------- # data mdata_spat <- op_spat %>% dplyr::select(Subject, Training_group, Trial, Rewarded, Choice_correct, Choice_number) %>% mutate(Rewarded=as.factor(Rewarded)) %>% mutate(Rewarded.lag=lag(Rewarded)) %>% droplevels() mdata_spat$Rewarded.lag[mdata_spat$Trial==1] <- "Y" # scale covariate mdata_spat$Trial.z <- scale(mdata_spat$Trial) # dummy-code and center factor variable mdata_spat$Rewarded.lag.Y <- (as.numeric(mdata_spat$Rewarded.lag)-1)-mean(as.numeric(mdata_spat$Rewarded.lag)-1) mdata_spat$Rewarded.Y <- (as.numeric(mdata_spat$Rewarded)-1)-mean(as.numeric(mdata_spat$Rewarded)-1) str(mdata_spat) # 320 obs # full model full_spat0 <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z |Subject), data = mdata_spat, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_spat0)$varcor # remove correlations between random intercept and random slopes full_spat <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z || Subject), data = mdata_spat, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Rewarded + Trial.z, data = mdata_spat) car::vif(xx) max(car::vif(xx)) # 1.006823 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_spat$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_spat[mdata_spat$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_spat@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_spat)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_spat <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_spat,3) # orig min max # (Intercept) -1.617 -2.167 -1.251 # Training_group5 -1.645 -2.111 -0.979 # RewardedY 0.581 0.306 0.719 # Trial.z 0.671 0.554 0.831 ggplot(mod.stab.res_full_spat, aes(x=orig, y = row.names(mod.stab.res_full_spat)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## not super stable, but ok # inference: .......................................................... null_spat <- glmer(Choice_correct ~ #Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z || Subject), data = mdata_spat, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_spat, null_spat, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_spat 6 244.37 266.79 -116.19 232.37 # full_spat 7 243.83 269.99 -114.92 229.83 2.5448 1 0.1107 ## no significant effect of group on correct choice drop1(full_spat, test="Chisq") # npar AIC LRT Pr(Chi) # 243.83 # Training_group 1 244.37 2.5448 0.11066 # Rewarded 1 243.77 1.9399 0.16368 # Trial.z 1 247.61 5.7798 0.01621 * round(summary(full_spat)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) -1.617 0.762 # Training_group5 -1.645 1.011 # RewardedY 0.581 0.386 # Trial.z 0.671 0.260 emmeans(full_spat, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.2061 0.1212 Inf 0.0573 0.526 # 5 0.0477 0.0337 Inf 0.0116 0.177 ## both groups are not above the chance level (0.16667) test(emmeans(full_spat, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.2061 0.1210 Inf 0.167 0.352 0.7246 # 5 0.0477 0.0337 Inf 0.167 -1.866 0.0621 ## both groups perform around chance in choosing the target bowl ## there is a trend for group 5 to choose below chance # model plot ................................................................... mdata_spat_plot <- as.data.frame(effects::effect("Training_group", full_spat)) ggplot(mdata_spat_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "OSI Test", x="Group", y="Probability of choosing the target bowl (2 or 5)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ranef.data_spat <- as.data.frame(ggeffects::ggpredict(full_spat, terms = c("Trial.z[all]", "Training_group", "Subject"), type="random")) names(ranef.data_spat) <- c("Trial.z", "predicted", "std.error", "conf.low", "conf.high", "Training_group", "Subject") ranef.data_spat <- ranef.data_spat %>% filter(Training_group=="2" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="2"]) | Training_group=="5" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="5"])) %>% mutate(Subject=factor(Subject, levels = c(sort(unique(op_spat$Subject[op_spat$Training_group=="2"])), sort(unique(op_spat$Subject[op_spat$Training_group=="5"])))), unz.trial = round((Trial.z * attr(mdata_spat$Trial.z, 'scaled:scale') + attr(mdata_spat$Trial.z, 'scaled:center'))),0) %>% filter(!(Subject=="Bibi"&unz.trial>=12)) %>% filter(!(Subject=="Zampano"&unz.trial>19)) ranef.data_spat %>% ggplot(aes(x = unz.trial, y = predicted))+ geom_hline(yintercept = 1/6, lty=2) + # chance level geom_line(aes(colour = Training_group), size=1) + # average probability geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Training_group), alpha=.2, show.legend = F) + geom_point(data=mdata_spat, aes(x = Trial, y = Choice_correct), inherit.aes = F, col="grey50") + facet_wrap(~Subject) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + scale_fill_manual(values=c("mediumorchid4","deepskyblue3")) + labs(title="Ordinal-Spatial Incongruence Test", x="Trials", y="Probability of choosing the target bowl (2 or 5)", col="Group") + theme_bw(base_size = 15) ## spatial 2/5 - first ----------------------------------------------- # data mdata_spat_f <- op_spat %>% dplyr::select(Subject, Training_group, Trial, Choice_correct, Choice_number) %>% filter(Trial<4) %>% droplevels() # scale covariate mdata_spat_f$Trial.z <- scale(mdata_spat_f$Trial) str(mdata_spat_f) # 48 obs # full model full_spat_f <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_spat_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_spat_f)$varcor # no need to remove correlations between random intercept and random slopes # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Trial.z, data = mdata_spat_f) car::vif(xx) max(car::vif(xx)) # 1 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_spat_f$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_spat_f[mdata_spat_f$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_spat_f@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(red.data, full.mod.stab, res, subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_spat_f)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_spat_f <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_spat_f,3) # orig min max # (Intercept) -8.178 -12.843 -2.698 # Training_group5 -0.824 -4.519 0.365 # Trial.z -0.238 -3.099 0.000 ggplot(mod.stab.res_full_spat_f, aes(x=orig, y = row.names(mod.stab.res_full_spat_f)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## not super stable, but to be expected with this little number of observations # inference: .......................................................... null_spat_f <- glmer(Choice_correct ~ #Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_spat_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_spat_f, null_spat_f, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_spat_f 5 35.563 44.919 -12.781 25.563 # full_spat_f 6 35.397 46.625 -11.699 23.397 2.1656 1 0.1411 ## no significant effect of group on correct choice drop1(full_spat_f, test="Chisq") # npar AIC LRT Pr(Chi) # 35.397 # Training_group 1 35.563 2.16564 0.1411 # Trial.z 1 33.413 0.01561 0.9006 round(summary(full_spat_f)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) -8.178 3.230 # Training_group5 -0.824 3.081 # Trial.z -0.238 1.931 emmeans(full_spat_f, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.000281 0.000906 Inf 5.00e-07 0.1362 # 5 0.000123 0.000418 Inf 1.58e-07 0.0875 ## in both groups the upper CI is below the chance level (0.16667) test(emmeans(full_spat_f, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.000281 0.000906 Inf 0.167 -2.034 0.0420 # 5 0.000123 0.000418 Inf 0.167 -2.177 0.0295 ## both groups perform significantly below chance in choosing the target bowl # model plot ................................................................... mdata_spat_f_plot <- as.data.frame(effects::effect("Training_group", full_spat_f)) ggplot(mdata_spat_f_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "OSI Test - First 3 Trials", x="Group", y="Probability of choosing the target bowl (2 or 5)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ## spatial 3/4 -------------------------------------------------- # data mdata_spat_34 <- op_spat %>% dplyr::select(Subject, Training_group, Trial, Rewarded, Choice_number) %>% mutate(Rewarded=as.factor(Rewarded)) %>% mutate(Rewarded.lag=lag(Rewarded)) %>% droplevels() mdata_spat_34$Rewarded.lag[mdata_spat$Trial==1] <- "Y" mdata_spat_34$Choice_correct <- ifelse(mdata_spat_34$Training_group=="2" & mdata_spat_34$Choice_number=="3" | mdata_spat_34$Training_group=="5" & mdata_spat_34$Choice_number=="4", 1, 0) # scale covariate mdata_spat_34$Trial.z <- scale(mdata_spat_34$Trial) # dummy-code and center factor variable mdata_spat_34$Rewarded.Y <- (as.numeric(mdata_spat_34$Rewarded)-1)-mean(as.numeric(mdata_spat_34$Rewarded)-1) str(mdata_spat_34) # 320 obs # full model full_spat_34 <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z |Subject), data = mdata_spat_34, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_spat_34)$varcor # no need to remove correlations between random intercept and random slopes # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Rewarded + Trial.z, data = mdata_spat_34) car::vif(xx) max(car::vif(xx)) # 1.006823 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_spat_34$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_spat_34[mdata_spat_34$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_spat_34@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_spat_34)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_spat_34 <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_spat_34,3) # orig min max # (Intercept) -0.386 -1.013 -0.062 # Training_group5 0.777 0.458 1.634 # RewardedY -0.400 -0.542 -0.280 # Trial.z -0.644 -0.729 -0.570 ggplot(mod.stab.res_full_spat_34, aes(x=orig, y = row.names(mod.stab.res_full_spat_34)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## group not very stable, but ok # inference: .......................................................... null_spat_34 <- glmer(Choice_correct ~ #Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z |Subject), data = mdata_spat_34, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_spat_34, null_spat_34, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_spat_34 9 381.87 415.49 -181.93 363.87 # full_spat_34 10 382.45 419.82 -181.23 362.45 1.4131 1 0.2346 ## no significant effect of group on correct choice drop1(full_spat_34, test="Chisq") # npar AIC LRT Pr(Chi) # 382.45 # Training_group 1 381.87 1.4131 0.2345502 # Rewarded 1 382.03 1.5743 0.2095894 # Trial.z 1 391.57 11.1217 0.0008533 *** round(summary(full_spat_34)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) -0.386 0.444 # Training_group5 0.777 0.619 # RewardedY -0.400 0.315 # Trial.z -0.644 0.165 emmeans(full_spat_34, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.362 0.106 Inf 0.188 0.583 # 5 0.553 0.107 Inf 0.347 0.742 ## both groups perform above the chance level (0.16667), see lower CI test(emmeans(full_spat_34, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.362 0.106 Inf 0.167 2.277 0.0228 # 5 0.553 0.107 Inf 0.167 4.227 <.0001 # model plot ................................................................. mdata_spat_34_plot <- as.data.frame(effects::effect("Training_group", full_spat_34)) ggplot(mdata_spat_34_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "OSI Test", x="Group", y="Probability of choosing the\nspatially correct bowl (3 or 4)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ranef.data_spat_34 <- as.data.frame(ggeffects::ggpredict(full_spat_34, terms = c("Trial.z[all]", "Training_group", "Subject"), type="random")) names(ranef.data_spat_34) <- c("Trial.z", "predicted", "std.error", "conf.low", "conf.high", "Training_group", "Subject") ranef.data_spat_34 <- ranef.data_spat_34 %>% filter(Training_group=="2" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="2"]) | Training_group=="5" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="5"])) %>% mutate(Subject=factor(Subject, levels = c(sort(unique(op_spat$Subject[op_spat$Training_group=="2"])), sort(unique(op_spat$Subject[op_spat$Training_group=="5"])))), unz.trial = round((Trial.z * attr(mdata_spat$Trial.z, 'scaled:scale') + attr(mdata_spat$Trial.z, 'scaled:center'))),0) %>% filter(!(Subject=="Bibi"&unz.trial>=12)) %>% filter(!(Subject=="Zampano"&unz.trial>19)) ranef.data_spat_34 %>% ggplot(aes(x = unz.trial, y = predicted))+ geom_hline(yintercept = 1/6, lty=2) + # chance level geom_line(aes(colour = Training_group), size=1) + # average probability geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Training_group), alpha=.2, show.legend = F) + geom_point(data=mdata_spat_34, aes(x = Trial, y = Choice_correct), inherit.aes = F, col="grey50") + facet_wrap(~Subject) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + scale_fill_manual(values=c("mediumorchid4","deepskyblue3")) + labs(title="Ordinal-Spatial Incongruence Test", x="Trials", y="Probability of choosing the spatially correct bowl (3 or 4)", col="Group") + theme_bw(base_size = 15) ## spatial 3/4 - first ----------------------------------------------- # data mdata_spat_34_f <- op_spat %>% dplyr::select(Subject, Training_group, Trial, Choice_number) %>% filter(Trial<4) %>% droplevels() mdata_spat_34_f$Choice_correct <- ifelse(mdata_spat_34_f$Training_group=="2" & mdata_spat_34_f$Choice_number=="3" | mdata_spat_34_f$Training_group=="5" & mdata_spat_34_f$Choice_number=="4", 1, 0) # scale covariate mdata_spat_34_f$Trial.z <- scale(mdata_spat_34_f$Trial) str(mdata_spat_34_f) # 48 obs # full model full_spat_34_f <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_spat_34_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_spat_34_f)$varcor # remove correlations between random intercept and random slopes, cite Matuschek 2017 full_spat_34_f_woc <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z ||Subject), data = mdata_spat_34_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) logLik(full_spat_34_f) # 'log Lik.' -26.53254 (df=6) logLik(full_spat_34_f_woc) # 'log Lik.' -26.92881 (df=5) ## log likelihood not very different (can use reduced model) ## woc model super unstable though so we keep the first model # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Trial.z, data = mdata_spat_34_f) car::vif(xx) max(car::vif(xx)) # 1 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_spat_34_f$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_spat_34_f[mdata_spat_34_f$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_spat_34_f@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_spat_34_f)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_spat_34_f <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_spat_34_f,3) # orig min max # (Intercept) 0.948 0.555 1.516 # Training_group5 0.288 -0.092 0.774 # Trial.z 1.416 1.211 1.904 ggplot(mod.stab.res_full_spat_34_f, aes(x=orig, y = row.names(mod.stab.res_full_spat_34_f)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## not super stable, but ok # inference: .......................................................... null_spat_34_f <- glmer(Choice_correct ~ #Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_spat_34_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) as.data.frame(anova(full_spat_34_f, null_spat_34_f, test="Chisq")) # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_spat_34_f 5 63.20983 72.56583 -26.60491 53.20983 NA NA NA # full_spat_34_f 6 65.06507 76.29228 -26.53254 53.06507 0.1447574 1 0.7035965 ## no significant effect of group on correct choice drop1(full_spat_34_f, test="Chisq") # npar AIC LRT Pr(Chi) # 65.065 # Training_group 1 63.210 0.1448 0.70360 # Trial.z 1 69.144 6.0794 0.01368 * round(summary(full_spat_34_f)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) 0.948 0.703 # Training_group5 0.288 0.759 # Trial.z 1.416 0.706 emmeans(full_spat_34_f, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.721 0.142 Inf 0.394 0.911 # 5 0.775 0.117 Inf 0.481 0.927 ## in both groups the lower CI is above the chance level (0.16667) test(emmeans(full_spat_34_f, pairwise ~ "Training_group", type="response"), null=boot::logit(1/6)) # Training_group prob SE df null z.ratio p.value # 2 0.721 0.142 Inf 0.167 3.637 0.0003 # 5 0.775 0.117 Inf 0.167 4.254 <.0001 ## both groups perform significantly above chance in choosing the spatial bowl # model plot ................................................................. mdata_spat_34_f_plot <- as.data.frame(effects::effect("Training_group", full_spat_34_f)) ggplot(mdata_spat_34_f_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/6, linetype = 2) + labs(title = "OSI Test - First 3 Trials", x="Group", y="Probability of choosing the\nspatially correct bowl (3 or 4)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ## fronto-parallel ----------------------------------------------------------- # data mdata_fp <- op_FP %>% dplyr::select(Subject, Training_group, Trial, Rewarded, Choice_correct, Choice_number) %>% mutate(Rewarded=as.factor(Rewarded)) %>% mutate(Rewarded.lag=lag(Rewarded)) %>% droplevels() mdata_fp$Rewarded.lag[mdata_fp$Trial==1] <- "Y" mdata_fp$Choice_number[mdata_fp$Choice_number==""] <- NA mdata_fp$Choice_correct[is.na(mdata_fp$Choice_number)] <- NA mdata_fp <- droplevels(mdata_fp) mdata_fp$search_bias <- ifelse(mdata_fp$Training_group=="2" & mdata_fp$Choice_number=="2L" | mdata_fp$Training_group=="5" & mdata_fp$Choice_number=="2R" , "left", ifelse(mdata_fp$Training_group=="5" & mdata_fp$Choice_number=="2L" | mdata_fp$Training_group=="2" & mdata_fp$Choice_number=="2R" , "right", "no" )) table(mdata_fp$search_bias, mdata_fp$Training_group) table(mdata_fp$Choice_correct, mdata_fp$Training_group) # scale covariate mdata_fp$Trial.z <- scale(mdata_fp$Trial) # dummy-code and center factor variable mdata_fp$Rewarded.Y <- (as.numeric(mdata_fp$Rewarded)-1)-mean(as.numeric(mdata_fp$Rewarded)-1) str(mdata_fp) # 320 obs # full model full_fp0 <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z |Subject), data = mdata_fp, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_fp0)$varcor # remove correlations between random intercept and random slopes full_fp <- glmer(Choice_correct ~ Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z ||Subject), data = mdata_fp, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Rewarded.lag + Trial.z, data = mdata_fp) car::vif(xx) max(car::vif(xx)) # 1.003082 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_fp$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_fp[mdata_fp$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_fp@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_fp)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_fp <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_fp,3) # orig min max # (Intercept) -0.284 -0.766 0.000 # Training_group5 0.472 0.179 0.942 # RewardedY 0.190 0.090 0.341 # Trial.z 0.626 0.488 0.728 ggplot(mod.stab.res_full_fp, aes(x=orig, y = row.names(mod.stab.res_full_fp)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## fine # inference: .......................................................... null_fp <- glmer(Choice_correct ~ #Training_group + Rewarded + Trial.z + (1 + Rewarded.Y + Trial.z ||Subject), data = mdata_fp, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_fp, null_fp, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_fp 6 371.18 393.37 -179.59 359.18 # full_fp 7 372.69 398.57 -179.35 358.69 0.4916 1 0.4832 ## no significant effect of group on correct choice drop1(full_fp, test="Chisq") # npar AIC LRT Pr(Chi) # 372.69 # Training_group 1 371.18 0.4916 0.483224 # Rewarded 1 371.16 0.4649 0.495354 # Trial.z 1 377.45 6.7593 0.009326 ** round(summary(full_fp)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) -0.284 0.521 # Training_group5 0.472 0.667 # RewardedY 0.190 0.276 # Trial.z 0.626 0.226 emmeans(full_fp, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.443 0.125 Inf 0.228 0.683 # 5 0.561 0.107 Inf 0.353 0.749 ## group 2 is not above the chance level (0.33333) in choosing 2L or 2R test(emmeans(full_fp, pairwise ~ "Training_group", type="response"), null=boot::logit(1/3)) # Training_group prob SE df null z.ratio p.value # 2 0.443 0.125 Inf 0.333 0.918 0.3588 # 5 0.561 0.107 Inf 0.333 2.165 0.0304 ## group 5 performs above chance level in choosing 2L or 2R # model plot ................................................................... mdata_fp_plot <- as.data.frame(effects::effect("Training_group", full_fp)) ggplot(mdata_fp_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/3, linetype = 2) + labs(title = "Fronto-Parallel Test", x="Group", y="Probability of choosing one of the \ntwo target bowls (2L-5R or 5L-2R)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ranef.data_fp <- as.data.frame(ggeffects::ggpredict(full_fp, terms = c("Trial.z[all]", "Training_group", "Subject"), type="random")) names(ranef.data_fp) <- c("Trial.z", "predicted", "std.error", "conf.low", "conf.high", "Training_group", "Subject") ranef.data_fp <- ranef.data_fp %>% filter(Training_group=="2" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="2"]) | Training_group=="5" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="5"])) %>% mutate(Subject=factor(Subject, levels = c(sort(unique(op_spat$Subject[op_spat$Training_group=="2"])), sort(unique(op_spat$Subject[op_spat$Training_group=="5"])))), unz.trial = round((Trial.z * attr(mdata_fp$Trial.z, 'scaled:scale') + attr(mdata_fp$Trial.z, 'scaled:center'))),0) %>% filter(!(Subject=="Bibi"&unz.trial>6)) %>% filter(!(Subject=="Raya"&unz.trial>=18)) %>% filter(!(Subject=="Zampano"&unz.trial>15)) ranef.data_fp %>% ggplot(aes(x = unz.trial, y = predicted))+ geom_hline(yintercept = 1/3, lty=2) + # chance level geom_line(aes(colour = Training_group), size=1) + # average probability # geom_line(aes(x = unz.trial, y = conf.low, colour = Training_group), size=1, lty=2) + # lower CI # geom_line(aes(x = unz.trial, y = conf.high, colour = Training_group), size=1, lty=2) + # upper CI geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Training_group), alpha=.2, show.legend = F) + geom_point(data=mdata_fp, aes(x = Trial, y = Choice_correct), inherit.aes = F, col="grey50") + facet_wrap(~Subject) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + scale_fill_manual(values=c("mediumorchid4","deepskyblue3")) + labs(title="Fronto-Parallel Test", x="Trials", y="Probability of choosing one of the two target bowls (2L-5R or 5L-2R)", col="Group") + theme_bw(base_size = 15) ## fronto-parallel - first ----------------------------------------------- # data mdata_fp_f <- op_FP %>% dplyr::select(Subject, Training_group, Trial, Choice_correct, Choice_number) %>% filter(Trial<4) %>% droplevels() mdata_fp_f$Choice_number[mdata_fp_f$Choice_number==""] <- NA mdata_fp_f$Choice_correct[is.na(mdata_fp_f$Choice_number)] <- NA mdata_fp_f <- droplevels(mdata_fp_f) # scale covariate mdata_fp_f$Trial.z <- scale(mdata_fp_f$Trial) str(mdata_fp_f) # 48 obs # full model full_fp_f0 <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z |Subject), data = mdata_fp_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_fp_f0)$varcor # remove correlations between random intercept and random slopes, cite Matuschek 2017 full_fp_f <- glmer(Choice_correct ~ Training_group + Trial.z + (1 + Trial.z ||Subject), data = mdata_fp_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_correct ~ Training_group + Trial.z, data = mdata_fp_f) car::vif(xx) max(car::vif(xx)) # 1 ## no issue with collinearity (should be around 1) ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_fp_f$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_fp_f[mdata_fp_f$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_fp_f@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } # mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(red.data, full.mod.stab, res, subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_fp_f)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_fp_f <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_fp_f,3) # orig min max # (Intercept) -1.138 -1.775 -0.886 # Training_group5 0.002 -0.624 0.444 # Trial.z -1.069 -1.620 -0.935 ggplot(mod.stab.res_full_fp_f, aes(x=orig, y = row.names(mod.stab.res_full_fp_f)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## ok # inference: .......................................................... null_fp_f <- glmer(Choice_correct ~ #Training_group + Trial.z + (1 + Trial.z ||Subject), data = mdata_fp_f, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_fp_f, null_fp_f, test="Chisq") # npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) # null_fp_f 4 57.935 65.420 -24.968 49.935 # full_fp_f 5 59.935 69.291 -24.968 49.935 0 1 0.9996 ## no significant effect of group on correct choice drop1(full_fp_f, test="Chisq") # npar AIC LRT Pr(Chi) # 59.935 # Training_group 1 57.935 0.000 0.99957 # Trial.z 1 64.554 6.619 0.01009 * round(summary(full_fp_f)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) -1.138 0.560 # Training_group5 0.002 0.739 # Trial.z -1.069 0.491 emmeans(full_fp_f, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.243 0.103 Inf 0.0965 0.490 # 5 0.243 0.102 Inf 0.0980 0.487 ## in both groups the lower CI is not above the chance level (0.3333) test(emmeans(full_fp_f, pairwise ~ "Training_group", type="response"), null=boot::logit(1/3)) # Training_group prob SE df null z.ratio p.value # 2 0.243 0.103 Inf 0.333 -0.794 0.4271 # 5 0.243 0.102 Inf 0.333 -0.802 0.4224 ## both groups perform around chance in choosing the target bowls (2L & 2R) in their first 3 trials # model plot .................................................................. mdata_fp_f_plot <- as.data.frame(effects::effect("Training_group", full_fp_f)) ggplot(mdata_fp_f_plot, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/3, linetype = 2) + labs(title = "Fronto-Parallel Test - First 3 Trials", x="Group", y="Probability of choosing one of the \ntwo target bowls (2L-5R or 5L-2R)") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) ## fronto-parallel - bowl bias ----------------------------------------- # data mdata_fp_bias <- mdata_fp %>% filter(Choice_correct==1) %>% droplevels() %>% mutate(search_bias=factor(search_bias, levels = c("right", "left"))) %>% mutate(search_bias_bin=as.numeric(as.factor(search_bias))-1, # right=0, left=1 Choice_number_bin=as.numeric(Choice_number)-1) # 2R=1, 2L=0 # get walking side per subject ws_sag <- op_sag %>% group_by(Subject) %>% count(Side, .drop = F) # check how many NA ws_sag[is.na(ws_sag$Side),] # only Barbarossa has NA # check what NA means for Barbarossa ws_sag[ws_sag$Subject=="Barbarossa",] # only left is missing # exchange NA for L ws_sag$Side[is.na(ws_sag$Side)] <- "L" # include missing rows ws_sag_exp <- expand.grid(Subject=unique(ws_sag$Subject), Side=unique(ws_sag$Side)) ws_sag_exp <- merge(ws_sag_exp, ws_sag, all.x = T) ws_sag_exp$n[is.na(ws_sag_exp$n)] <- 0 # remove straight trials and calculate proportions ws_sag_Rprop <- ws_sag_exp %>% filter(Side!="straight") %>% group_by(Subject) %>% mutate(walk.side_prop=prop.table(n)) %>% filter(Side=="R") # merge with bias data set mdata_fp_wbias <- merge(mdata_fp_bias, ws_sag_Rprop[, c("Subject", "walk.side_prop")]) # z-transform new covariate mdata_fp_wbias$walk.side_prop.z <- scale(mdata_fp_wbias$walk.side_prop) str(mdata_fp_wbias) # 151 obs ## full model full_fp_wbias0 <- glmer(Choice_number_bin ~ Training_group + walk.side_prop.z + Trial.z + (1 + Trial.z |Subject), data = mdata_fp_wbias, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) summary(full_fp_wbias0)$varcor # remove correlations between random intercept and random slopes full_fp_wbias <- glmer(Choice_number_bin ~ Training_group + walk.side_prop.z + Trial.z + (1 + Trial.z ||Subject), data = mdata_fp_wbias, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # assumptions ......................................................... ## collinearity: xx <- lm(Choice_number_bin ~ Training_group + Trial.z + walk.side_prop.z, data = mdata_fp_wbias) car::vif(xx) max(car::vif(xx)) # 1.1687 ## no issue with collinearity (should be around 1) ## group and walk side are somewhat related ## model stability: library(foreach) library(doParallel) # How many cores does your CPU have n_cores <- detectCores() cluster <- makeCluster(n_cores - 1) # Register cluster registerDoParallel(cluster) subjects <- unique(as.character(mdata_fp_wbias$Subject)) i=1 mod.stab.res_sub <- foreach(i = 1:length(subjects), .combine=cbind) %dopar% { # remove 1 level from data set red.data <- mdata_fp_wbias[mdata_fp_wbias$Subject!=subjects[i],] red.data <- droplevels(red.data) # run model with reduced data set full.mod.stab <- lme4::glmer(formula = full_fp_wbias@call$formula, data = red.data, family = binomial(link="logit"), control = lme4::glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) # get coefficients res <- as.data.frame(summary(full.mod.stab)$coefficients[,1]) names(res) <- subjects[i] res # must be run because the last object in the loop will be used to cbind to results object } mod.stab.res_sub # stop the cluster! stopCluster(cl = cluster) # remove helper objects rm(subjects, i) # get final result mod.stab.res_orig <- as.data.frame(summary(full_fp_wbias)$coefficients[,1]) names(mod.stab.res_orig) <- "orig" mod.stab.res <- cbind(mod.stab.res_orig, mod.stab.res_sub) mod.stab.res$min <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = min, na.rm=TRUE) mod.stab.res$max <- apply(X = mod.stab.res[,-1], MARGIN = 1, FUN = max, na.rm=TRUE) mod.stab.res_full_fp_wbias <- mod.stab.res %>% dplyr::select(orig, min, max) round(mod.stab.res_full_fp_wbias,3) # orig min max # (Intercept) 3.454 2.501 5.492 # Training_group5 -0.565 -2.527 0.672 # walk.side_prop.z 1.845 1.243 2.264 # Trial.z 0.406 0.217 0.598 ggplot(mod.stab.res_full_fp_wbias, aes(x=orig, y = row.names(mod.stab.res_full_fp_wbias)))+ geom_point() + geom_errorbarh(aes(xmin = min, xmax = max), height = .2) + geom_vline(xintercept = 0, linetype = 2) + theme_bw() ## ok # inference: .......................................................... null_fp_wbias <- glmer(Choice_number_bin ~ # Training_group + walk.side_prop.z + Trial.z + (1 + Trial.z || Subject), data = mdata_fp_wbias, family = binomial(link="logit"), control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))) anova(full_fp_wbias, null_fp_wbias, test="Chisq") # npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) # null_fp_wbias 5 103.42 118.50 -46.708 93.416 # full_fp_wbias 6 105.32 123.43 -46.662 93.324 0.092 1 0.7616 ## no significant effect of group on bowl choice drop1(full_fp_wbias, test="Chisq") # npar AIC LRT Pr(Chi) # 105.32 # Training_group 1 103.42 0.0920 0.76163 # walk.side_prop.z 1 108.27 4.9421 0.02621 * # Trial.z 1 104.90 1.5792 0.20887 ## effect of walking side in sag test on choice in fp test round(summary(full_fp_wbias)$coefficients[,1:2],3) # Estimate Std. Error # (Intercept) 3.454 1.539 # Training_group5 -0.565 1.846 # walk.side_prop.z 1.845 1.061 # Trial.z 0.406 0.325 ## pigs that had a right walking bias, also chose 2R more likely emmeans(full_fp_wbias, pairwise~"Training_group", type="response")$emmeans # Training_group prob SE df asymp.LCL asymp.UCL # 2 0.980 0.0292 Inf 0.722 0.999 # 5 0.868 0.1169 Inf 0.470 0.980 test(emmeans(full_fp_wbias, pairwise ~ "Training_group", type="response"), null=boot::logit(1/2)) # Training_group prob SE df null z.ratio p.value # 2 0.971 0.0429 Inf 0.5 2.286 0.0223 # 5 0.951 0.0638 Inf 0.5 2.175 0.0297 ## both groups significantly prefer 2R over 2L if walkside is included in model ### model plots .............................................................. # choice per group mdata_fp_wbias_plot1 <- as.data.frame(effects::effect("Training_group", full_fp_wbias)) ggplot(mdata_fp_wbias_plot1, aes(x = Training_group, y = fit, color = Training_group)) + geom_point(size=3, show.legend = F) + geom_errorbar(aes(ymin = lower, ymax = upper), width = .2, show.legend = F) + geom_hline(yintercept = 1/2, linetype = 2) + labs(title = "Fronto-Parallel Test - Bowl Bias", x="Group", y="Probability of choosing 2R over 2L") + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + ylim(0,1) + theme_bw(base_size = 15) # choice over walkingside mdata_fp_wbias_plot2 <- as.data.frame(ggeffects::ggpredict(full_fp_wbias, "walk.side_prop.z")) mdata_fp_wbias_plot2$walk.side_prop <- round(mdata_fp_wbias_plot2$x * attr(mdata_fp_wbias$walk.side_prop.z, 'scaled:scale') + attr(mdata_fp_wbias$walk.side_prop.z, 'scaled:center'),5) # raw data points mdata_fp_wbias_plot2_raw <- mdata_fp_wbias %>% group_by(Subject, walk.side_prop) %>% summarise(Choice_number_bin.prop = mean(Choice_number_bin)) ggplot(mdata_fp_wbias_plot2, aes(x = walk.side_prop, y = predicted)) + geom_point(data = mdata_fp_wbias_plot2_raw, mapping = aes(x = walk.side_prop, y = Choice_number_bin.prop), alpha = .2, size=3) + geom_smooth(se = F) + geom_smooth(aes(x = walk.side_prop, y = conf.low), se = F, linetype = 2) + geom_smooth(aes(x = walk.side_prop, y = conf.high), se = F, linetype = 2) + geom_hline(yintercept = 1/2, linetype = 2) + labs(title = "Transfer Test - Bowl Bias", x="Proportion of Trials with Right Walking Side\nin the Sagittal Test", y="Probability of choosing 2R over 2L") + ylim(0,1) + xlim(0,1)+ theme_bw(base_size = 15) # per subject ranef.data_fp_wbias <- as.data.frame(ggeffects::ggpredict(full_fp_wbias, terms = c("Trial.z[all]", "Training_group", "Subject"), type="random")) names(ranef.data_fp_wbias) <- c("Trial.z", "predicted", "std.error", "conf.low", "conf.high", "Training_group", "Subject") ranef.data_fp_wbias <- ranef.data_fp_wbias %>% filter(Training_group=="2" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="2"]) | Training_group=="5" & Subject%in%unique(op_spat$Subject[op_spat$Training_group=="5"])) %>% mutate(Subject=factor(Subject, levels = c(sort(unique(op_spat$Subject[op_spat$Training_group=="2"])), sort(unique(op_spat$Subject[op_spat$Training_group=="5"])))), unz.trial = round((Trial.z * attr(mdata_fp$Trial.z, 'scaled:scale') + attr(mdata_fp$Trial.z, 'scaled:center')), 0)) ranef.data_fp_wbias <- ranef.data_fp_wbias %>% filter(!(Subject=="Bibi")) %>% filter(!(Subject=="Blossom"&unz.trial<3), !(Subject=="Blossom"&unz.trial>19)) %>% filter(!(Subject=="Raya"&unz.trial>15)) %>% filter(!(Subject=="Zazou"&unz.trial<7)) %>% filter(!(Subject=="Zoltan"&unz.trial<14)) %>% filter(!(Subject=="Bessy"&unz.trial>19)) %>% filter(!(Subject=="Bijou"&unz.trial<6)) %>% filter(!(Subject=="Radomir"&unz.trial>18)) %>% filter(!(Subject=="Ronon"&unz.trial<4)) %>% filter(!(Subject=="Zampano"&unz.trial>6)) ranef.data_fp_wbias %>% ggplot(aes(x = unz.trial, y = predicted))+ geom_hline(yintercept = 1/2, lty=2) + # chance level geom_line(aes(colour = Training_group), size=1) + # average probability geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Training_group), alpha=.2, show.legend = F) + geom_point(data=mdata_fp_wbias, aes(x = Trial, y = ifelse(Choice_number=="2R",1,0)), inherit.aes = F, col="grey50") + facet_wrap(~Subject) + scale_color_manual(values=c("mediumorchid4","deepskyblue3")) + scale_fill_manual(values=c("mediumorchid4","deepskyblue3")) + labs(title="Fronto-Parallel Test", x="Trials", y="Probability of choosing 2R over 2L", col="Group") + theme_bw(base_size = 15)