####ANALYSIS ABSOLUTELY RELATIVE CHICKS#### ####Load Packages#### library(ggplot2) library(dplyr) library(multimode) library(lme4) library(glmmTMB) library(emmeans) ####Read data and define variable classes#### discr <- read.csv("Data Absolutely Relative Chicks.csv", sep=";", dec=",", tryLogical=F) #last one ensures that sex does not become true or false str(discr) discr <- discr[discr$Subject<125,] discr$Subject <- as.factor(discr$Subject) discr$Homogenous_Heterogenous <- as.factor(discr$Homogenous_Heterogenous) discr$Sex <- as.factor(discr$Sex) discr$Training <- as.factor(discr$Training) discr$Positive_Training_Numerosity <- as.factor(discr$Positive_Training_Numerosity) discr$Positive_Higher_Lower <- as.factor(discr$Positive_Higher_Lower) discr$Test <- as.factor(discr$Test) discr$Side_Positive <- as.factor(discr$Side_Positive) discr$Zone_Left <- as.numeric(discr$Zone_Left) discr$Zone_Middle <- as.numeric(discr$Zone_Middle) discr$Zone_Right <- as.numeric(discr$Zone_Right) discr$Zone_Correct <- as.numeric(discr$Zone_Correct) discr$Zone_Incorrect <- as.numeric(discr$Zone_Incorrect) discr$Group <- as.factor(discr$Group) ####Create a data set per individual#### discr_subj <- discr %>% group_by(Subject) %>% summarise( Zone_Correct_mean = mean(Zone_Correct), Zone_Correct = sum(Zone_Correct), Zone_Incorrect = sum(Zone_Incorrect), Zone_Middle = sum(Zone_Middle), Sex, Training, Positive_Training_Numerosity, Test, Side_Positive, Group, Preference_Index=mean(Preference_Index)) #take only every 6th row discr_subj <- discr_subj[seq(1, nrow(discr_subj ),6),] ####calculate sample size per experiment#### summary(discr_subj) ####Calculate cumulative preference index per subject#### discr_subj$Preference_Index_cum <- discr_subj$Zone_Correct/(discr_subj$Zone_Correct+discr_subj$Zone_Incorrect) summary(discr_subj) Group1 <- subset(discr_subj, Group==1) Group2 <- subset(discr_subj, Group==2) Group3 <- subset(discr_subj, Group==3) Group4 <- subset(discr_subj, Group==4) #control group summary(discr$Subject) summary(discr$Preference_Index) ###Testing against chance level (group-level strategy preference)#### #Preference index per minute m1 <- glmmTMB(Preference_Index ~ 1 + Group + (1|Subject), ordbeta(link="logit"), data=discr) summary(m1) #Intercept n.s. #Confidence intervals per group per minute emmeans(m1, ~Group, type="response") #Confidence intervals per group and subject emmeans(m2, ~Group, type="response") #preference index per subject m2 <- glmmTMB(Preference_Index ~ 1 + Group, ordbeta(link="logit"), data=discr_subj) summary(m2, specs="Group") #intercept n.s. #confidence intervals preference index per minute: emmip(m1,~Group, type="response", engine="ggplot", CIs=T, CIarg=list(lwd=1.5, alpha=0.2), xlab= "Group", ylab= "Preference Index", linearg=list(col="white", alpha=0), dotarg=list(size=4, col="turquoise")) + theme_classic() + scale_y_continuous(limits=c(0,1)) + theme(axis.title=element_text(size=15, face="bold"), axis.text=element_text(size=10)) + geom_hline(aes(yintercept=0.5), linetype="dashed", size=1) #confidence intervals preference index per subject emmip(m2,~Group, type="response", engine="ggplot", CIs=T, CIarg=list(lwd=1.5, alpha=0.2), xlab= "Group", ylab= "Preference Index", linearg=list(col="white", alpha=0), dotarg=list(size=4, col="turquoise")) + theme_classic() + scale_y_continuous(limits=c(0,1)) + theme(axis.title=element_text(size=15, face="bold"), axis.text=element_text(size=10)) + geom_hline(aes(yintercept=0.5), linetype="dashed", size=1) ###Differences between experiments 1 and 2#### wilcox.test(Group1$Preference_Index, Group2$Preference_Index) #W = 865.5, p-value = 0.5315 ####Investigating the individual preference distributions#### #main density plot ggplot(data=discr_subj, aes( x=Preference_Index_cum, group=Group, fill=Group, alpha=0.2)) + geom_density(data=discr_subj, stat="density") + theme_classic() + xlab("Preference Index per Subject") + ylab("Density") + scale_alpha(guide=F) #last part removes alpha from the legend ####Testing for bimodality (against unimodality)#### #Ameijeiras-Alonso et al. (2019) excess mass test #Group1 modetest(Group1$Preference_Index_cum) #Excess mass = 0.14304, p-value = 0.038 locmodes(Group1$Preference_Index_cum, mod0=2) #Estimated location #Modes: 0.3733176 0.9575151 #Antimode: 0.682744 modetest(Group2$Preference_Index_cum) #Excess mass = 0.16082, p-value = 0.008 locmodes(Group2$Preference_Index_cum, mod0=2) #Estimated location #Modes: 0.367636 0.9302031 #Antimode: 0.6725359 modetest(Group3$Preference_Index_cum) #Excess mass = 0.22668, p-value = 0.002 locmodes(Group3$Preference_Index_cum, mod0=2) #Estimated value of the density #Estimated location #Modes: 0.02331283 0.5553202 #Antimode: 0.2792899 modetest(Group4$Preference_Index_cum) #Excess mass = 0.174, p-value = 0.122 #Comparing the distributions of the individual preference indices between experiments 1 and 2 ks.test(Group1$Preference_Index_cum, Group2$Preference_Index_cum) #D = 0.175, p-value = 0.5607 ####Effect of the position (left or right) of the reinforced stimulus on the strategy preference#### #full model m3 <- glmmTMB(Preference_Index ~ Side_Positive*Group + (1|Subject), ordbeta(link="logit"), data=discr) summary(m3) #null model m3_null <- glmmTMB(Preference_Index ~ 1 + (1|Subject), ordbeta(link="logit"), data=discr) anova(m3, m3_null) #p=0.2664 #full model m4 <- glmmTMB(Preference_Index ~ Side_Positive*Group, ordbeta(link="logit"), data=discr_subj) summary(m4) #null model m4_null <- glmmTMB(Preference_Index ~ 1, ordbeta(link="logit"), data=discr_subj) anova(m4, m4_null) #p=0.569 #plots: #all ggplot(data=discr, aes(x=Group, y=Preference_Index)) + geom_boxplot(aes(fill=Side_Positive)) + theme_classic() + ggtitle("PI (per min and subject) ~ Group*Side") #individual ggplot(data=discr_subj, aes(x=Group, y=Preference_Index)) + geom_boxplot(aes(fill=Side_Positive)) + theme_classic() + ggtitle("PI per subject ~ Group*Side")