##### Preliminaries ############################################################

# Install and load required packages. We'll be using ordinal-logistic regression
# via the clm() function: a specific type of generalised linear model.
install.packages(c("ordinal", "sure"))
library(ordinal)
library(sure)

#load data Exp 1
library(readxl)

#Adjust script to set the correct working directory
DataExp1 <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp1.xlsx")
View(DataExp1)


# Convert numerics in to factors.
DataExp1$Block <- as.factor(DataExp1$Block)
DataExp1$Tool <- as.factor(DataExp1$Tool)

############################  Model 1  ##################################
## Test if jays (i) adjust their preference according to the apparatus 
#(Wide vs Narrow tube); (ii) adjust their preference across SessionBlocks
##(e.g. acquire a preference for the small stone)

# Fit ordinal-logit model to data from Experiment 1 (stones as tool).
#Include terms to account for Block and Apparatus.
model1 <- clm(Tool ~ Block  + Tube,
              data = subset(DataExp1, Experiment == "1A"))
summary(model1)

# Find that:
#    - proportion of selection of the 3 stones is sig dif only in 1 case Block1 vs Block4; but not sig
#      for all other comparisons (i.e. Bl 1 vs bl 2; Bl 1 vs Bl 3)
#    - The Apparatus has no significant effect on stone choice.

############################  Model 2  ##################################
## Test if jays transfer their preference from stones to novel objects 
# We compare the Block 4 (last 20 trials with stones (Experiment 1A)) with Block 5
#(the first 20 trials with novel objects (Experiment 1B ))
#NB 1 bird (Chinook) must be excluded because she did not participate to Exp 1B

# Fit ordinal-logit model to data from birds other than Chinook, from Blocks 4
# and 5. Include terms to account for Experiment, Subject, and Apparatus.
model2 <- clm(Tool ~ Tube + Experiment,
              data = subset(DataExp1, Subject != "CHINOOK" &
                              Block %in% c(4, 5)))
summary(model2)
# Find that:
#    - proportion of selection of the 3 stones is not different between Block 4 and block 5
#    - The Apparatus has no significant effect on choice.

############################  Model 3  ##################################

## Test if jays preference changes across Blocks (5 and 6) in Exp1B (e.g. if
#the percentage of selection of the small tool continue to increase in Exp2.

# Fit ordinal-logit model to data from Blocks 5 and 6. Include terms to account
# for Block, Subject, and Apparatus.
model3 <- clm(Tool ~ Block  + Tube,
              data = subset(DataExp1, Block %in% c(5, 6)))
summary(model3)


# Find that:

#    - Block has a significant effect on choice; the birds appear to choose
#      smaller objects in Block 6.
#    - The Apparatus has no significant effect on choice.

##### Model Diagnostics ########################################################

# Diagnostics for ordinal regression models are complicated. However, for each
# of the clm() models, we can use autoplot.clm() from the sure package to
# produce some useful diagnostic plots. For example, for the model fitted above,
# we would enter
autoplot.clm(model1, what = "qq")
autoplot.clm(model2, what = "qq")
autoplot.clm(model3, what = "qq")
# and
autoplot.clm(model1, what = "fitted")
autoplot.clm(model2, what = "fitted")
autoplot.clm(model3, what = "fitted")
#############################################################################
#############################################################################

##############################     EXP 2     ##########################


## GLM
########## run a GLM to test whether the two stones are selected in different proportions 

#load data
library(readxl)
DataExp2 <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp2.xlsx")
View(DataExp2)


#convert numerics into factors
DataExp2$Subject       <- as.factor(DataExp2$Subject)
DataExp2$Tube        <- as.factor(DataExp2$Tube)
DataExp2$Block        <- as.factor(DataExp2$Block)
DataExp2$Tool       <- as.factor(DataExp2$Tool)

##GLM
model4 <- glm(Tool ~ Tube  + Block,
               family = binomial(link = "logit"), data = DataExp2)

summary(model4)

## Find that, preference for Long stone:
## 1.  is not dependent on the conditions
## 2.  is stable across the 2 sessions

##########################################################
##################### Tool Manipulation


########### GLM 
library(readxl)
DataExp2Rot <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp2.xlsx", 
                       sheet = "ToolRotation")
View(DataExp2Rot)


#convert numerics into factors
DataExp2Rot$Rotation2Lev        <- as.factor(DataExp2Rot$Rotation2Lev)
DataExp2Rot$Subject       <- as.factor(DataExp2Rot$Subject)
DataExp2Rot$Tube        <- as.factor(DataExp2Rot$Tube)
DataExp2Rot$Block        <- as.factor(DataExp2Rot$Block)

##GLM
model5 <- glm(Rotation2Lev ~ Tube  + Block,
               family = binomial(link = "logit"), data = DataExp2Rot)

summary(model5)

# Find that
# 1. jays's manipulation is sig different between condition (Wide vs Narrow Tube)
# 2. there is no effect of block

########################################################
########################################################
########################################################

## BOXPLOT Exp1
library(readxl)
library(ggplot2)

DataExp1Plot <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp1Plot.xlsx")
View(DataExp1Plot)

# convert numerics into factor & Re-order 'Tool' (so that the order is: S->M->L).
DataExp1Plot$Tool <- factor(DataExp1Plot$Tool, levels = c("Small", "Medium", "Large"))
DataExp1Plot$Block <- factor(DataExp1Plot$Block)

# grouped boxplot
ggplot(DataExp1Plot, aes(x=Block, y=PropSel, fill=Tool)) + 
  geom_boxplot() + theme_classic() + scale_fill_manual(values=c("#E69F00",  "#56B4E9", "#999999"))

#######################################################################
#######################################################################
#### Individual Plots Exp 1a-b

library(readxl)
library(ggplot2)
library(dplyr)

#load data
DataExp1 <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp1.xlsx")
View(DataExp1) 

# create subset data for each individual
Chinook = subset(DataExp1, Subject == "CHINOOK")
Homer = subset(DataExp1, Subject == "HOMER")
Jaylo = subset(DataExp1, Subject == "JAYLO")
Poe = subset(DataExp1, Subject == "POE")
Stuka = subset(DataExp1, Subject == "STUKA")


### Chinook
ggplot( Chinook, aes(x=Trial, y=Tool, colour=Tube)) +
  geom_line( color="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 80, 10)) +
  scale_y_discrete(limits=c("1","2","3"), labels=c("Small", "Medium", "Large"),expand = c(0.05,0.05))

### HOMER
ggplot( Homer, aes(x=Trial, y=Tool, colour=Tube)) +
  geom_line( color="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 120, 10)) +
  scale_y_discrete(limits=c("1","2","3"), labels=c("Small", "Medium", "Large"), expand = c(0.05,0.05))

### JAYLO
ggplot( Jaylo, aes(x=Trial, y=Tool, colour=Tube)) +
  geom_line( color="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 120, 10)) +
  scale_y_discrete(limits=c("1","2","3"), labels=c("Small", "Medium", "Large"), expand = c(0.05,0.05))

### POE
ggplot( Poe, aes(x=Trial, y=Tool, colour=Tube)) +
  geom_line( color="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 120, 10)) +
  scale_y_discrete(limits=c("1","2","3"), labels=c("Small", "Medium", "Large"), expand = c(0.05,0.05))

### STUKA
ggplot( Stuka, aes(x=Trial, y=Tool, colour=Tube)) +
  geom_line( color="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 120, 10)) +
  scale_y_discrete(limits=c("1","2","3"), labels=c("Small", "Medium", "Large"), expand = c(0.05,0.05))

######################################################################################################
######################################################################################################
######################################################################################################

#### Individual Plots Exp 2
library(readxl)
library(ggplot2)
library(dplyr)

#Load Data
DataExp2Plots <- read_excel("Desktop/ToolUseJaysPaper/Zenodo/DataExp2Plots.xlsx")
View(DataExp2Plots)

#create subset data for each individual
Homer2 = subset(DataExp2Plots, Subject == "HOMER")
Jaylo2 = subset(DataExp2Plots, Subject == "JAYLO")
Poe2 = subset(DataExp2Plots, Subject == "POE")
Stuka2 = subset(DataExp2Plots, Subject == "STUKA")


### HOMER2
ggplot( Homer2, aes(x=Trial, y=TooL_Rot_Code, colour=Tube)) +
  geom_line(colour="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 20, 5)) +
  scale_y_discrete(limits=c("Round","Long_N","Long_E","Long_I"), expand = c(0.05,0.05))


### Jaylo2
ggplot( Jaylo2, aes(x=Trial, y=TooL_Rot_Code, colour=Tube)) +
  geom_line(colour="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 20, 5)) +
  scale_y_discrete(limits=c("Round","Long_N","Long_E","Long_I"), expand = c(0.05,0.05))

### POE2
ggplot( Poe2, aes(x=Trial, y=TooL_Rot_Code, colour=Tube)) +
  geom_line(colour="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 20, 5)) +
  scale_y_discrete(limits=c("Round","Long_N","Long_E","Long_I"), expand = c(0.05,0.05))

### Stuka2
ggplot( Stuka2, aes(x=Trial, y=TooL_Rot_Code, colour=Tube)) +
  geom_line(colour="grey") +
  geom_point(aes(colour = Tube), size=2) + scale_color_manual(values=c("#56B4E9", "#E69F00"))+
  theme_classic()+  scale_x_continuous(breaks=seq(0, 20, 5)) +
  scale_y_discrete(limits=c("Round","Long_N","Long_E","Long_I"), expand = c(0.05,0.05))
######################################################################