################################### LIBRARY ##############################################

library(brms)
library(carData)
library(car)
library(bayestestR)
library(ggplot2)
library(viridis)
library(parallel)
library(gridExtra)
library(survminer)
library(ggpubr)
library(cowplot)
library(harrypotter)
detectCores()

####################################IMPORT DATASET#########################################

xdata=read.table(file="S1_Data", header=T, sep=";")
xdata
head(xdata)
table(xdata$task)
summary(xdata)

##############################  1. DEVELOPMENT #############################################################
ddata=subset(xdata, action=="insert")
table(ddata$grasp.type) #FH=29, FHT=243, GC=1188
table(ddata$subject) #70 individuals

#CHECK DATA#

table(ddata$subject, ddata$sex) 
hist(ddata$age) 
table(ddata$age, ddata$subject)
table(ddata$subject, ddata$group)
table(ddata$grasp.type, ddata$group)
table(ddata$age, ddata$grasp.type)
table(ddata$grasp.type)
table(ddata$group)
table(ddata$context)

## Z-TRANSFORM COVARIATE ##

z.age.d=as.vector(scale(ddata$age))
ddata$response<-relevel(as.factor(ddata$grasp.type), ref="FH")
ddata=data.frame(ddata, z.age.d)
rm(z.age.d)

##CHECK FOR COLLINEARITY##

cresponse.d = as.numeric (ddata$response)
cmodel.d=lm(cresponse.d ~ z.age.d + sex + group, data=ddata)
vif(cmodel.d)		

####MODEL####
ddata$response<-relevel(as.factor(ddata$grasp.type), ref="FH")

test.development<- brm (response ~ z.age.d + I(z.age.d^2) + sex + group + 
                          (1|subject) + (1|session) + (1+z.age.d+I(z.age.d^2)|context), 
                        data=ddata, family="categorical",cores=12, warmup = 1000, 
                        iter = 2000, chains = 4, control=list(adapt_delta = 0.999, max_treedepth=12, stepsize=0.01))

pp_check(test.development, ndraws=100)
summary(test.development)
ci(test.development, ci=0.89)
ci(test.development, ci=0.95)

####PLOT####
names(ddata)

zdata<-ddata[,c(7,26,20)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "age", "grasp", "Num")

a=aggregate(zdata$Num, by=list(zdata$age, zdata$subject), sum)
b=aggregate(zdata$Num, by=list(zdata$age, zdata$subject, zdata$grasp), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("age","subject","grasp","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total 

legend_labels <- c("Digits Grip", "Full Hand Grip", "Full Hand Thumb Grip")

d.plot<-conditional_effects(test.development, categorical=TRUE, effects="z.age.d")

#graph
((1:55)-mean(ddata$age))/sd(ddata$age) #for breaks and labels

plot.test.dev<-plot(d.plot, plot = FALSE)[[1]] +
  labs(y="Probability to use each hand grip\n")+
  scale_y_continuous(limits=c(0,1))+
  geom_point(data=b, mapping=aes(x=age, y=proba, size=number, colour=grasp, shape=grasp) ,alpha=0.65, 
             inherit.aes = FALSE)+
  geom_line(linewidth=3)+
  scale_size_continuous(name="number", range=c(3, 20))+
  scale_x_continuous(name ="\nAge (years)", 
                     limits=NULL,
                     breaks=c(-1.50352561, -1.19653433,-0.81279522,-0.42905612, -0.04531707, 0.33842209, 0.72216119, 1.10590030,1.48963940, 1.87337851, 2.25711761,2.64085672),
                     labels=c("1", "5", "10", "15", "20", "25","30","35", "40", "45", "50", "55"))+
  scale_fill_manual(values=c("black","orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("FH"="Full Hand Grip", "FHT"="Full Hand Thumb Grip", "DG"="Digits Grip"))+                #Change couleur CI
  scale_color_manual(values=c("black", "orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("FH"="Full Hand Grip", "FHT"="Full Hand Thumb Grip", "DG"="Digits Grip"))+ 
  theme_classic()+
  guides(shape = guide_legend(ncol=1, override.aes = list(size = 10)), fill = guide_legend(ncol = 1), size = guide_legend(ncol = 1))+
  scale_shape_discrete(limits=c("FH", "FHT", "DG"),labels=c("FH"="Full Hand Grip", "FHT"="Full Hand Thumb Grip", "DG"="Digits Grip"))+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="top",
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))


plot.test.dev

################################## ESTIMATE OF GRASP SEEN ON PLOT ####################################
amc <- conditional_effects(test.development, categorical = TRUE)
amc_data <- amc[[1]]
amc_data

names(amc_data)
estimate.age <- amc_data[,c(11,12,13)]
colnames(estimate.age)= c("age", "grasp", "estimate")

write.csv2(estimate.age, "estimate grasp per age.csv", row.names = FALSE)

##############################  2. ACCURACY #############################################################
adata= subset(ddata, attempts!="Unclear")
table(adata$attempts)
adata$attempts=as.numeric(adata$attempts)
nrow(adata)

#CHECK DATA#
table(adata$attempts)
hist(adata$age)
table(adata$attempts, adata$sex)
table(adata$group, adata$attempts)

# Z-TRANSFORM COVARIATE #
z.age.a=as.vector(scale(adata$age))
adata=data.frame(adata,z.age.a)
rm(z.age.a)
adata$grasp.type<-relevel(as.factor(adata$grasp.type), ref="FH")

##CHECK FOR COLLINEARITY##

cmodel.a=lm(attempts ~ grasp.type + z.age.a + sex + group, data=adata)
vif(cmodel.a)		

####MODEL####
adata$grasp.type<-relevel(as.factor(adata$grasp.type), ref="FH")

test.accuracy<- brm (attempts ~ grasp.type + z.age.a + I(z.age.a^2) + sex + group +
                       (1|subject) + (1|session) + (1 + z.age.a +I(z.age.a^2) |context), 
                     data=adata, family="poisson", cores=12, warmup = 1000, 
                     iter = 2000, chains = 4, control=list(adapt_delta = 0.999, max_treedepth=12, stepsize=0.01))

pp_check(test.accuracy, ndraws=100)
summary(test.accuracy)
ci(test.accuracy, ci=0.89)
ci(test.accuracy, ci=0.95)

#RELEVEL MODEL#
adata$grasp.type<-relevel(as.factor(adata$grasp.type), ref="FHT")
test.accuracy.r<- brm (attempts ~ grasp.type + z.age.a + I(z.age.a^2) + sex + group +
                         (1|subject) + (1|session) + (1 + z.age.a +I(z.age.a^2) |context), 
                       data=adata, family="poisson",cores=12, warmup = 1000, 
                       iter = 2000, chains = 4, control=list(adapt_delta = 0.999, max_treedepth=12, stepsize=0.01))

pp_check(test.accuracy.r, ndraws=100)
summary(test.accuracy.r)
ci(test.accuracy.r, ci=0.89)
ci(test.accuracy.r, ci=0.95)

####PLOT####
names(adata)

zdata<-ddata[,c(7,24,20)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "attempts", "grasp", "Num")

a=aggregate(zdata$Num, by=list(zdata$grasp), sum)
b=aggregate(zdata$Num, by=list(zdata$grasp, zdata$attempts), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("grasp","attempts","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total
b$attempts=as.numeric(b$attempts)

a.plot <- (conditional_effects(test.accuracy, categorical=FALSE, effects="grasp.type", re_formula=NA,
                               robust=TRUE, method="fitted",
                               spaghetti=FALSE, surface=FALSE, resolution=100))

plot.test.acc <- plot(a.plot, plot= FALSE) [[1]] +                  
  theme(plot.title=element_text(hjust=0.5))+  
  xlab("\nHand grip type") + ylab("Number of attempts\n") +
  geom_boxplot(aes(color=grasp.type))+
  geom_point(data=b, mapping=aes(x=grasp, y=attempts,  size=number, colour=grasp) ,alpha=0.4,
             inherit.aes = FALSE)+
  scale_size_continuous(name="number", range=c(3, 20))+
  scale_fill_manual (values=c("black","orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("FH"="Full hand grip", "FHT"="Full hand thumb grip", "DG"="Digits grip"))+
  scale_color_manual(values=c("black","orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("FH"="Full hand grip", "FHT"="Full hand thumb grip", "DG"="Digits grip"))+
  scale_x_discrete (limits=c("FH", "FHT", "DG"),labels=c("FH"="Full hand grip", "FHT"="Full hand thumb grip", "DG"="Digits grip"))+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="top",
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))

plot.test.acc

##############################  3. ACTION: INSERT VS POUND #############################################################
############# 3.1 ADULT INDIVIDUALS ##########################
data3.1=subset(xdata, action != "lever")
data3.2=subset(data3.1, action != "prod")
data3.3=subset(data3.2, action != "screw")
data3.4=subset(data3.3, action != "stir")
data3.5=subset(data3.4, pounding.ind=="Yes") ## Only individuals that do pounding action in the dataset
tdata=subset(data3.5, age>15)

tdata$action<-relevel(as.factor(tdata$action), ref="insert")
tdata$response<-relevel(as.factor(tdata$grasp.type), ref="DG")

#CHECK DATA#
table(tdata$action)
table(tdata$grasp.type)
hist(tdata$age)
table(tdata$grasp.type, tdata$action)

#CHECK COLLINEARITY
cresponse.t = as.numeric (tdata$response)
cmodel.t=lm(cresponse.t ~ action + sex + group, data=tdata)
vif(cmodel.t)	

##### MODEL #####
tdata$action<-relevel(as.factor(tdata$action), ref="insert")
tdata$response<-relevel(as.factor(tdata$grasp.type), ref="DG")

bprior <- prior(normal(0,1), class= b)

test.action.pound.insert.adult<- brm (response ~ action + sex + group +
                                        (1|subject) + (1|session) + (1|context), 
                                      data=tdata, family="categorical", prior=bprior, warmup = 1000, 
                                      iter = 2000, chains = 4, cores=12, control=list(adapt_delta = 0.999, max_treedepth=15))

pp_check(test.action.pound.insert.adult, ndraws=100)
summary(test.action.pound.insert.adult)
ci(test.action.pound.insert.adult, ci=0.89)
ci(test.action.pound.insert.adult, ci=0.95)

####PLOT####
names(tdata)

zdata<-tdata[,c(6,15,14)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "action", "grasp",  "Num")
zdata$grasp<- factor(zdata$grasp, c("FH", "FHT", "DG"))

a=aggregate(zdata$Num, by=list(zdata$task, zdata$subject), sum)
b=aggregate(zdata$Num, by=list(zdata$task, zdata$subject, zdata$grasp), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("action","subject","grasp","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total 

t.plot.pound <- (conditional_effects(test.action.pound.insert.adult, categorical=TRUE, effects="action", re_formula=NA,
                                     robust=TRUE, method="fitted",
                                     spaghetti=FALSE, surface=FALSE, resolution=100))

plot_new<-t.plot.pound
plot_new$`task:cats__`$effect2__ <- factor(plot_new$`task:cats__`$effect2__, c("FH", "FHT", "DG"))

plot.test.action.pound.adult <- plot(plot_new, plot= FALSE) [[1]] +                  
  theme(plot.title=element_text(hjust=0.5))+  
  xlab("\nAction") + ylab("Probability to use each hand grip\n") +
  geom_boxplot()+
  geom_point(data=b, mapping=aes(x=action, y=proba,  size=number, colour=grasp) ,alpha=0.2,
             inherit.aes = FALSE,position = position_dodge(width = .4))+
  scale_size_continuous(name="number", range=c(3, 20))+
  scale_fill_manual (values=c("black","orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("DG"="Digits grip","FH"="Full hand grip", "FHT"="Full hand thumb grip"))+
  scale_color_manual (values=c("black","orange", "#56B4E9"), limits=c("FH", "FHT", "DG"),labels=c("DG"="Digits grip","FH"="Full hand grip", "FHT"="Full hand thumb grip"))+
  scale_x_discrete (limits=c("insert", "pounding"),labels=c("insert"="Inserting","pound"="Pounding"))+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="top",
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))

plot.test.action.pound.adult

############# 3.2 NON-ADULT INDIVIDUALS ##########################
tdata2=subset(data3.5, age<15)

tdata2$action<-relevel(as.factor(tdata2$action), ref="insert")
tdata2$response<-relevel(as.factor(tdata2$grasp.type), ref="DG")

nrow(tdata2)

#CHECK DATA#
table(tdata2$action)
table(tdata2$grasp.type)
hist(tdata2$age)
table(tdata2$grasp.type, tdata2$action)

#CHECK COLLINEARITY
cresponse.t2 = as.numeric (tdata2$response)
cmodel.t2=lm(cresponse.t2 ~ action + sex + group, data=tdata2)
vif(cmodel.t2)	

##### MODEL #####
tdata2$action<-relevel(as.factor(tdata2$action), ref="insert")
tdata2$response<-relevel(as.factor(tdata2$grasp.type), ref="DG")

bprior <- prior(normal(0,1), class= b)

test.action.pound.insert.NONadult<- brm (response ~ action + sex + group +
                                           (1|subject) + (1|session) + (1|context), 
                                         data=tdata2, family="categorical", prior=bprior, warmup = 1000, 
                                         iter = 2000, chains = 4, cores=12, control=list(adapt_delta = 0.999, max_treedepth=15))

pp_check(test.action.pound.insert.NONadult, ndraws=100)
summary(test.action.pound.insert.NONadult)
ci(test.action.pound.insert.NONadult, ci=0.89)
ci(test.action.pound.insert.NONadult, ci=0.95)

##############################  4. LARVAE EXTRACTION - ACTION SUCCESS #############################################################
data4.1=subset(xdata, action != "insert")
data4.2=subset(data4.1, larvae_success != "Unclear")
ldata=subset(data4.2, context =="larvae_extraction")

#CHECK DATA#
table(ldata$action)
table(ldata$context)
hist(ldata$age)
table(ldata$age.class)
table(ldata$subject, ldata$sex)
table(ldata$sex)
table(ldata$levering)
table(ldata$group)
table(ldata$larvae_success, ldata$action)

#CHECK COLLINEARITY#
ldata$response<-relevel(as.factor(ldata$larvae_success), ref="Whole")
cresponse.l = as.numeric (ldata$response)
cmodel.t=lm(cresponse.l ~ age.class + levering + sex, data=ldata)
vif(cmodel.t)	

##MODEL##
ldata$response<-relevel(as.factor(ldata$larvae_success), ref="None")

bprior <- prior(normal(0,1), class= b)

test.larvae.success <- brm (response ~ age.class + levering + sex +
                              (1|subject) + (1|session), 
                            data=ldata, family="categorical", prior=bprior, warmup = 1000,  
                            iter = 2000, chains = 4, cores=3, control=list(adapt_delta = 0.999, max_treedepth=15, stepsize=0.01))

pp_check(test.larvae.success, ndraws=100)
summary(test.larvae.success)
ci(test.larvae.success, ci=0.89)
ci(test.larvae.success, ci=0.95)

####PLOT####
names(ldata)

zdata<-ldata[,c(6,20,19)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "outcome", "levering", "Num")
zdata$outcome<- factor(zdata$outcome, c("None", "Partial", "Whole"))

a=aggregate(zdata$Num, by=list(zdata$outcome, zdata$subject), sum)
b=aggregate(zdata$Num, by=list(zdata$outcome, zdata$subject, zdata$levering), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("outcome","subject","levering","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total 

plot.levering.success<-conditional_effects(test.larvae.success, categorical=TRUE, effects="levering")

plot_new<-plot.levering.success
plot_new$`levering:cats__`$effect2__ <- factor(plot_new$`levering:cats__`$effect2__, c("None", "Partial", "Whole"))


plot.test.levering.success <- plot(plot_new, plot= FALSE) [[1]] +                  
  theme(plot.title=element_text(hjust=0.5))+  
  xlab("\nAction") + ylab("Probability to retrieve \ndifferent amount of larvae\n") +
  geom_boxplot()+
  geom_point(data=b, mapping=aes(x=levering, y=proba,  size=number, colour=outcome) ,alpha=0.4,
             inherit.aes = FALSE, position = position_dodge(width = .4))+
  scale_size_continuous(name="number", range=c(3, 20))+
  scale_fill_hp_d (option="LunaLovegood", limits=c("None", "Partial", "Whole"),labels=c("None"="None", "Partial"="Partial", "Whole"="Whole"))+
  scale_color_hp_d(option="LunaLovegood", limits=c("None", "Partial", "Whole"),labels=c("None"="None", "Partial"="Partial", "Whole"="Whole"))+
  scale_x_discrete (limits=c("yes", "no"),labels=c("no"="Other", "yes"="Levering"))+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="right",
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))

plot.test.levering.success

##############################  5. AGE-CLASS USE OF LEVERING ACTION #############################################################
#CHECK COLLINEARITY#
ldata$response<-relevel(as.factor(ldata$levering), ref="no")
cresponse.l = as.numeric (ldata$response)
cmodel.t=lm(cresponse.l ~ age.class + sex, data=ldata)
vif(cmodel.t)	

##MODEL##
ldata$response<-relevel(as.factor(ldata$levering), ref="no")

test.levering.ageclass <- brm (response ~ age.class + sex +
                                 (1|subject) + (1|session), 
                               data=ldata, family="bernoulli", prior=bprior, warmup = 1000,  
                               iter = 2000, chains = 4, cores=3, control=list(adapt_delta = 0.999, max_treedepth=15, stepsize=0.01))

pp_check(test.levering.ageclass, ndraws=100)
summary(test.levering.ageclass)
ci(test.levering.ageclass, ci=0.89)
ci(test.levering.ageclass, ci=0.95)

####PLOT####
names(ldata)

zdata<-ldata[,c(6,10,19)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "age.class", "levering", "Num")

a=aggregate(zdata$Num, by=list(zdata$age.class, zdata$subject), sum)
b=aggregate(zdata$Num, by=list(zdata$age.class, zdata$subject, zdata$levering), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("age.class","subject","levering","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total 

plot.levering.age.class<-conditional_effects(test.levering.ageclass, categorical=FALSE, effects="age.class")


plot.test.levering.age.class <- plot(plot.levering.age.class, plot= FALSE) [[1]] +                  
  theme(plot.title=element_text(hjust=0.5))+  
  xlab("\nAge-class") + ylab("Probability to use\nthe levering action\n") +
  geom_boxplot(aes(color=age.class))+
  geom_point(data=b, mapping=aes(x=age.class, y=proba,  size=number, colour=age.class) ,alpha=0.4,
             inherit.aes = FALSE)+
  scale_size_continuous(name="number", range=c(3, 20))+
  scale_fill_hp_d (option="Gryffindor", limits=c("adult", "non-adult"),labels=c("adult"="Adult individuals", "non-adult"="Non-adult individuals"))+
  scale_color_hp_d(option="Gryffindor", limits=c("adult", "non-adult"),labels=c("adult"="Adult individuals", "non-adult"="Non-adult individuals"))+
  scale_x_discrete (limits=c("adult", "non-adult"),labels=c("adult"="Adult individuals\n\u2265 15 years old", "non-adult"="Non-adult individuals\n < 15 years old"))+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="right",
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))

plot.test.levering.age.class

##############################  6. LEVERING ACTION IN DIFFERENT CONTEXTS #############################################################
############# 6.1 LARVAE EXTRACTION CONTEXT ##########################
data6.1=subset(xdata, action != "insert")
data6.2=subset(data6.1, age<20)

ldata.l=subset(data6.2, context =="larvae_extraction")

#CHECK DATA#
table(ldata.l$action)
table(ldata.l$context)
hist(ldata.l$age)
table(ldata.l$subject, ldata$sex)
table(ldata.l$sex)
table(ldata.l$levering)
table(ldata.l$group)

#CHECK COLLINEARITY
z.age.l=as.vector(scale(ldata.l$age))
ldata.l=data.frame(ldata.l, z.age.l)
rm(z.age.l)

#CHECK COLLINEARITY#
ldata.l$levering=as.factor(ldata.l$levering)
ldata.l$levering=as.numeric(ldata.l$levering)
cmodel.t=lm(levering ~ z.age.l + sex, data=ldata.l)
vif(cmodel.t)	

##### MODEL #####
ldata.l$levering<-relevel(as.factor(ldata.l$levering), ref="no")

bprior <- prior(normal(0,1), class= b)

test.larvae.spoon <- brm (levering ~ z.age.l +
                            (1|subject) + (1|session), 
                          data=ldata.l, family="bernoulli", warmup = 1000,  
                          iter = 2000, chains = 4, cores=12, control=list(adapt_delta = 0.99, max_treedepth=15))

pp_check(test.larvae.spoon, ndraws=100)
summary(test.larvae.spoon)
ci(test.larvae.spoon, ci=0.87)
ci(test.larvae.spoon, ci=0.89)
ci(test.larvae.spoon, ci=0.95)

####PLOT####
names(ldata.l)

zdata<-ldata[,c(6,22,19)]
zdata$Num=1
head(zdata)
colnames(zdata)=c("subject", "age", "levering", "Num")

a=aggregate(zdata$Num, by=list(zdata$age, zdata$subject), sum)
b=aggregate(zdata$Num, by=list(zdata$age, zdata$subject, zdata$levering), sum)

b$total_context=NA

for (i in 1:nrow(b)){
  xx=which(a$Group.1 == b$Group.1[i] & a$Group.2==b$Group.2[i])
  if(length(xx)>0){
    b$total_context[i]=as.character(a$x[xx[1]])}
}

colnames(b)=c("age","subject","levering","number","total")
b$number=as.numeric(b$number)
b$total=as.numeric(b$total)
b$proba=b$number/b$total 

plot.larvae.spoon<-conditional_effects(test.larvae.spoon, categorical=FALSE, effects="z.age.l", prob=0)
((1:55)-mean(ldata$age))/sd(ldata$age) #for breaks and labels

plot.test.l.spoon<-plot(plot.larvae.spoon, plot = FALSE)[[1]] +
  labs(y="Probability to use\n the levering action\n")+
  scale_y_continuous(limits=c(0,1))+
  scale_x_continuous(name ="\nAge (years)", 
                     limits=c(-1.57638289, 1.36604523),
                     breaks=c(-1.57638289, -0.95692434, -0.18260115, 0.59172204, 1.36604523),
                     labels=c("1", "5", "10", "15", "20"))+
  theme_classic()+
  geom_line(linewidth=3, color="aquamarine4")+
  geom_point(data=b, mapping=aes(x=age, y=proba, size=number) ,alpha=0.4, color="aquamarine4",
             inherit.aes = FALSE)+
  theme(legend.title=element_blank(), legend.text = element_text(size=35), legend.position="right", 
        axis.title.x = element_text(size=40, face="bold"),
        axis.title.y = element_text(size=40, face="bold"),
        axis.text.x = element_text(size=40, face="bold"),
        axis.text.y = element_text(size=40, face="bold"))

plot.test.l.spoon

############# 6.2 NUT KERNEL EXTRACTION CONTEXT ##########################
ndata=subset(data6.2, context =="nut_extraction")

#CHECK DATA#
table(ndata$action)
table(ndata$context)
hist(ndata$age)
table(ndata$subject, ndata$sex)
table(ndata$sex)
table(ndata$levering)
table(ndata$group)

z.age.n=as.vector(scale(ndata$age))
ndata=data.frame(ndata, z.age.n)
rm(z.age.n)

#CHECK COLLINEARITY#
ndata$levering=as.factor(ndata$levering)
ndata$levering=as.numeric(ndata$levering)
cmodel.t=lm(levering ~ z.age.n+sex, data=ndata)
vif(cmodel.t)	

##### MODEL #####
ndata$levering<-relevel(as.factor(ndata$levering), ref="no")

bprior <- prior(normal(0,1), class= b)

test.nut.spoon <- brm (levering ~ z.age.n +
                         (1|subject) + (1|session), 
                       data=ndata, family="bernoulli", warmup = 1000,  
                       iter = 2000, chains = 4, cores=3, control=list(adapt_delta = 0.999, max_treedepth=15))

pp_check(test.nut.spoon, ndraws=100)
summary(test.nut.spoon)
ci(test.nut.spoon, ci=0.85)
ci(test.nut.spoon, ci=0.89)
ci(test.nut.spoon, ci=0.95)

####PLOT####
names(ndata)
zdata2<-ndata[,c(6,22,19)]
zdata2$Num=1
head(zdata2)
colnames(zdata2)=c("subject", "age", "spooning", "Num")

a2=aggregate(zdata2$Num, by=list(zdata2$age, zdata2$subject), sum)
b2=aggregate(zdata2$Num, by=list(zdata2$age, zdata2$subject, zdata2$spooning), sum)

b2$total_context=NA

for (i in 1:nrow(b2)){
  xx=which(a2$Group.1 == b2$Group.1[i] & a2$Group.2==b2$Group.2[i])
  if(length(xx)>0){
    b2$total_context[i]=as.character(a2$x[xx[1]])}
}

colnames(b2)=c("age","subject","spooning","number","total")
b2$number=as.numeric(b2$number)
b2$total=as.numeric(b2$total)
b2$proba=b2$number/b2$total 

plot.nut.spoon<-conditional_effects(test.nut.spoon, categorical=FALSE, effects="z.age.n")

((1:55)-mean(ndata$age))/sd(ndata$age) #for breaks and labels

plot.test.n.spoon<-plot(plot.nut.spoon, plot = FALSE)[[1]] +
  labs(y="Probability to use levering in nut kernel extraction\n")+
  scale_y_continuous(limits=c(0,1))+
  scale_x_continuous(name ="\nAge", 
                     limits=c(-1.96149633,1.38843641),
                     breaks=c(-1.96149633, -1.25624733,-0.37468609,0.50687516, 1.38843641),
                     labels=c("1", "5", "10", "15", "20"))+
  theme_classic()+
  geom_line(linewidth=1.5)+
  geom_point(data=b2, mapping=aes(x=age, y=proba, size=number) ,alpha=0.3,
             inherit.aes = FALSE)+
  theme(legend.title=element_blank(), legend.text = element_text(size=30), legend.position="top", 
        axis.title.x = element_text(size=35, face="bold"),
        axis.title.y = element_text(size=35, face="bold"),
        axis.text.x = element_text(size=35, face="bold"),
        axis.text.y = element_text(size=35, face="bold"))

plot.test.n.spoon


############# PLOT - MERGE NUT KERNEL AND LARVAE EXTRACTION PLOTS TOGETHER ##################
#Put nut kernel regression line on same scale as larvae
0-mean(ldata$age)/sd(ldata$age)
vraix1=(plot.test.n.spoon$data$effect1__*sd(ndata$age))+mean(ndata$age)
vraix2=(vraix1-mean(ldata$age))/sd(ldata$age)
y1=plot.nut.spoon$z.age.n$estimate__

#extract larvae CI
plot.larvae.spoon1<-conditional_effects(test.larvae.spoon, categorical=FALSE, effects="z.age.l")
ymin1=plot.larvae.spoon1$z.age.l$lower__
ymax1=plot.larvae.spoon1$z.age.l$upper__

#put datapoint of nut kernel on same scale as larvae
age.nut=(b2$age*sd(ndata$age))+mean(ndata$age)
b2=data.frame(b2, age.nut)
z.age.mod=(b2$age.nut-mean(ldata$age))/sd(ldata$age)
b2=data.frame(b2, z.age.mod)

#create a new plot name for larvae extraction
plot<- plot.test.l.spoon
#add nut kernel datapoints and regression line to it
plot<- plot + geom_point(data=b2, mapping=aes(x=z.age.mod, y=proba, size=number) , color="burlywood4", alpha=0.4, inherit.aes = FALSE,show.legend = F)+
  geom_ribbon(aes(ymin=plot.test.n.spoon$data$lower__, ymax=plot.test.n.spoon$data$upper__, x=vraix2), alpha=0.1, 
              color = "burlywood4", linetype = "solid", fill="burlywood")
#add CI of larvae to it
plot<- plot + geom_ribbon(aes(x=vraix2, ymin=y1, ymax=y1), alpha=0.1, fill = NA, 
                          color = "burlywood4", linetype = "solid", size=3)
#add CI of nut kernel to it
plot <- plot + geom_ribbon(aes(ymin=ymin1, ymax=ymax1), alpha=0.1, fill = "aquamarine3", 
                           color = "aquamarine4", linetype = "solid")
#uniformization of colors
colors <- c("Nut kernel extraction"="burlywood4", "Larvae extraction"="aquamarine4")
#make it prettier by having the datapoints bigger
plot <- plot + geom_line(aes(y=plot.test.l.spoon$data$estimate__, color = "Larvae extraction"), size = 1.5) +
  geom_point (aes(x=vraix2, y=y1, color= "Nut kernel extraction"), size=1.5)+
  scale_colour_manual(name="Context",values=colors, 
                      guide = guide_legend(override.aes=aes(fill=NA)))+
  scale_size_continuous(name="number", range=c(5, 12))

plot