###### Food preference study with herring gull chicks ######
##Corresponding author: Emma Inzani eli204@exeter.ac.uk


##Influence of diet treatment group, test day and individual identity on food preferences

### Please load lines 4-118 
# Preliminary data tests 120-159
# Graphs lines 161-278
# Statistics line 280-496
# Inter Observer comparison stats lines 498- 695
# Chick biometrics statistics and graphs lines 696- END

dframe1 <- read.csv("Rehab_Chicks Foodpreference.csv")
head(dframe1) #puts an x in front of all numbers in column names
dim(dframe1[1]) #145 obs
names(dframe1)
library(ggplot2)
library(lme4)
library(performance)
library(MASS)
library(nnet)
library(ordinal)

# Remove Non Study birds
dframe1 <- dframe1[dframe1$Treatment!="Nonstud",]
dim(dframe1[1]) #n=135 obs

summary(dframe1)
class(dframe1$first.food.eaten)


#Changing variable classes to factors
dframe1$Chick.ID <- as.factor(dframe1$Chick.ID)
dframe1$Test <- as.factor(dframe1$Test)
dframe1$first.food.eaten <- as.factor(dframe1$first.food.eaten)
dframe1$second.food <- as.factor(dframe1$second.food)
dframe1$First.choice.front.row<- as.factor(dframe1$First.choice.front.row)
dframe1$First.choice.right <- as.factor(dframe1$First.choice.right)
dframe1$Second.choice.front<- as.factor(dframe1$Second.choice.front)
dframe1$second.choice.right <- as.factor(dframe1$second.choice.right)
dframe1$familiar.diet <- as.factor(dframe1$familiar.diet)
dframe1$Food.type <- as.factor(dframe1$Food.type)
dframe1$high.protein <- as.factor(dframe1$high.protein)
dframe1$secondhigh.protein <- as.factor(dframe1$secondhigh.protein)



#Mapping First and second food choice variables to recognise NAs
levels(dframe1$first.food.eaten)[levels(dframe1$first.food.eaten)==""] <- NA
dframe1$first.food.eaten
levels(dframe1$second.food)[levels(dframe1$second.food)==""] <- NA
dframe1$second.food
levels(dframe1$First.choice.front.row)[levels(dframe1$First.choice.front.row)==""] <- NA
dframe1$First.choice.front.row
levels(dframe1$First.choice.right)[levels(dframe1$First.choice.right)==""] <- NA
dframe1$First.choice.right
levels(dframe1$Second.choice.front)[levels(dframe1$Second.choice.front)==""] <- NA
dframe1$Second.choice.front
levels(dframe1$second.choice.right)[levels(dframe1$second.choice.right)==""] <- NA
dframe1$second.choice.right


#Reorder Test day lvls
dframe1$Test <- factor(dframe1$Test , levels=c("day1" ,"day5", "day10", "day15", "day35"))

##Subsetting individual tests
dframe2 <- dframe1[dframe1$Test=="day5",]
dim(dframe2[1]) #n=27 chicks -  only 25 made food choices
dframe2$first.food.eaten
head(dframe2)
summary(dframe2)
dframe3 <- dframe1[dframe1$Test=="day10",]
dframe4 <- dframe1[dframe1$Test=="day15",]
dframe5 <- dframe1[dframe1$Test=="day35",]
dframe20 <- dframe1[dframe1$Test=="day1",]

##Removing NAs from first food choice variables
dframe6 <- dframe2[complete.cases(dframe2$first.food.eaten),] # 5 day
dframe7 <- dframe3[complete.cases(dframe3$first.food.eaten),] # 10 day
dframe8 <- dframe4[complete.cases(dframe4$first.food.eaten),] # 15 day
dframe9 <- dframe5[complete.cases(dframe5$first.food.eaten),] # 35 day
dframe24 <- dframe20[!is.na(dframe20$first.food.eaten),] # baseline


dim(dframe6[1]) # n= 25
dim(dframe7[1]) # n= 26
dim(dframe8[1]) # n= 24
dim(dframe9[1]) # n= 13
dim(dframe24[1])# n= 23

## Removing NAs from second food choice for each test trial
dframe26 <- dframe2[complete.cases(dframe2$second.food),] # 5 day
dframe27 <- dframe3[complete.cases(dframe3$second.food),] # 10 day
dframe28 <- dframe4[complete.cases(dframe4$second.food),] # 15 day
dframe29 <- dframe5[complete.cases(dframe5$second.food),] # 35 day

dim(dframe26[1]) # n= 24
dim(dframe27[1]) # n= 22
dim(dframe28[1]) # n= 21
dim(dframe29[1]) # n= 7


#Full dataframe complete.cases subsets
dframe100 <- dframe1[dframe1$Test!="day1",] #removing 24hr test that doesn't have full array
dframe10 <- dframe100[complete.cases(dframe100$first.food.eaten),] #complete cases without 24hr
dframe101 <- dframe1[complete.cases(dframe1$first.food.eaten),] #including 24hr
#complete cases second choice without 24hr
dframe110 <- dframe100[complete.cases(dframe100$second.food),] 

dim(dframe1[1]) # n= 135 all test trials conducted including day1 test trials (unsuccessful and successful trials)
dim(dframe101[1]) # n= 111 successful test trials only, including day1 test trials
dim(dframe100[1]) # n= 108 all test trials conducted without day1 test trials (unsuccessful and successful trials)

dim(dframe10[1]) # n=88, successful only test trials, without day1 test trials as not directly comparable
dim(dframe110[1]) # n= 74 successful only second food choices, without day1 test trials



### Checking First and Second choices bias front/back right/left food array

#All trials but 24hr Front of array
table1 = table(dframe10$First.choice.front.row) # NA removed
table1 
chisq.test(table1) # X-squared = 49.5, df = 1, p-value = 1.984e-12 front bias,
table1 = table(dframe10$First.choice.front.row,dframe10$first.food.eaten) 
table1 
chisq.test(table1) #X^2= 9.55 df=3 P= 0.023 choice was also important, choose fish
fisher.test(table1)

table2 = table(dframe110$Second.choice.front) 
table2
chisq.test(table2) # X-squared = 2.6486, df = 1, p-value = 0.1036
table2 = table(dframe110$Second.choice.front,dframe110$second.food) 
table2 #X-squared = 10.11, df = 3, p-value = 0.01765
fisher.test(table2)

#Right side including 24hr=dframe101
table1 = table(dframe101$First.choice.right) 
table1 #exactly 44:44 for dframe10, 52:59 for dframe101
chisq.test(table1) # X^2= 0  P=1 dframe10, X-squared = 0.44144, df = 1, p-value = 0.5064
table1 = table(dframe101$First.choice.right,dframe101$first.food.eaten) 
table1 
chisq.test(table1) #X-squared = 1.1733, df = 3, p-value = 0.7594
fisher.test(table1)
# dframe 101 X-squared = 3.6091, df = 3, p-value = 0.3069

table2 = table(dframe10$second.choice.right) 
table2
chisq.test(table2) # X-squared = 0.21622, df = 1, p-value = 0.6419
table2 = table(dframe10$second.choice.right,dframe10$second.food) 
table2 #X^2= 3.08 P=0.38 no bias
chisq.test(table2) #X-squared = 4.7288, df = 3, p-value = 0.1928
fisher.test(table2)
###Only front of array bias also with preference for fish, no side bias###



####Visualizing data####
########################

####### Boxplot of latency to eat first food vs age/experience in captivity ($Test) ########
ggplot(dframe10, aes(x=Test, y=log(food.latency+1))) + 
  geom_boxplot(outlier.colour="black", outlier.shape=16,outlier.size=2, fill=c("#D55E00", "#E69F00", "#56B4E9", "#009E73"))+
  geom_jitter(position=position_jitter(0.2), color="gray55", shape=4, cex=2) +
  labs(x="Test day during captivity", y="Log latency to consume food (s)") +
  theme(legend.text=element_text(size=10),legend.title=element_text(size=0)) + theme_bw() + theme(strip.background =element_rect(fill="gray100"))+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))


###################### First food eaten graphs ######################

#### Test day and first foods eaten only ####
#############################################
##Figure 1
ggplot(dframe10, aes(x=Test, fill=first.food.eaten, facet = first.food.eaten))+
  geom_bar(position= position_dodge2(preserve = "single"), width=1)+
  scale_fill_manual(values = c("#663300", "#E69F00", "#56B4E9", "#009E73"), labels = c("Bread", "Cat food", "Fish", "Mussels")) +
  facet_wrap(~Treatment) +
  scale_y_continuous(name="Count", limits=c(0, 10), breaks =seq(0,10,2))+
  scale_x_discrete(name = "Test day",labels = c("Day5","Day10","Day15","Day35")) +
  labs(x="Food eaten first", y="Frequency of choice") +
  guides(fill=guide_legend(title="Food eaten first")) + theme_bw() + theme(text = element_text(size=25)) +theme(legend.text=element_text(size=20),legend.title=element_text(size=20)) + 
   theme(strip.background =element_rect(fill="gray85"))+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))



#### Individual chick first choices ####
######################################## 
##Figure 2
a<- ggplot(dframe10, aes(x=Chick.ID, fill=first.food.eaten)) + 
  geom_bar(position= "stack", width=0.7) +  
  labs(x="Chick IDs", y="Frequency of choice") +
  scale_fill_manual(values = c("#663300", "#E69F00", "#56B4E9", "#009E73"), labels = c("Bread", "Cat food", "Fish", "Mussels")) +
  theme_bw() + theme(text = element_text(size=17)) + theme(legend.position="none") +guides(fill="none") + theme(strip.background =element_rect(fill="gray100")) + guides(fill=guide_legend(title="First food eaten")) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
annotate(geom="text", x=4, y=4.5, label="A", color="black", size =10)
#most individuals have only 3 choices recorded out of 4 trials, 3 individuals have 2 choices
a

b<- ggplot(dframe110, aes(x=Chick.ID, fill=second.food)) + 
  geom_bar(position= "stack", width=0.7) +  
  labs(x="Chick IDs", y=) +
  scale_fill_manual(values = c("#663300", "#E69F00", "#56B4E9", "#009E73"), labels = c("Bread", "Cat food", "Fish", "Mussels")) +
  theme_bw() + theme(text = element_text(size=17)) + theme(strip.background =element_rect(fill="gray100")) + guides(fill=guide_legend(title="Food Choice")) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  annotate(geom="text", x=4, y=4.5, label="B", color="black", size =10)#
#most individuals have only 3 choices recorded out of 4 trials, 3 individuals have 2 choices
b

library(patchwork)
a+b
#### Amounts eaten of different foods by Treatment group ####
#############################################################
##Figure 3

#Need to alter dataframe shape
library(tidyr)
#make two columns with measure and second with values, all columns to be converted in list
datalong <- gather(dframe10, key="measure", value="value", c("amount.fish.eaten", "amount.mussel.eaten", "amount.cat.food.eaten","amount.bread.eaten"))
datalong
#make key variable and rename variables within
variable_names <- list(
  "amount.fish.eaten" = "Fish",
  "amount.mussel.eaten" = "Mussels",
  "amount.cat.food.eaten" = "Cat food",
  "amount.bread.eaten" = "Bread"
)

group_names <- levels(dframe10$Treatment)

variable_labeller <- function(variable,value){
  return(variable_names[value])
}

variable_labeller2 <- function(variable,value){
  if (variable=='measure') {
    return(variable_names[value])
  } else {
    return(group_names)
  }
}

#Changing name of Treatment groups
datalong$Treatment
mapping <- c("Marine" = "Marine", "Terrestrial" = "Terrestrial")
datalong$Treatment <- mapping[datalong$Treatment]

##Graph
ggplot(datalong, aes(x=Test, y=value, color=Treatment))+
  geom_boxplot(outlier.colour="black", outlier.shape=16,outlier.size=1.25)+
  geom_jitter(position=position_jitter(0.2), shape=4, cex=1.25) +
  facet_wrap(~measure, labeller = variable_labeller) +
  scale_color_manual(values = c(Marine = "#56B4E9", Terrestrial = "#E69F00" )) +
  labs(x="Test day", y="Amount eaten (g)") +
  scale_y_continuous(name="Amount eaten (g)", limits=c(0, 10), breaks =seq(0,10,2))+
  theme_bw() + theme(text = element_text(size=25)) +  theme(strip.background =element_rect(fill="gray"))+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))




################# Protein of foods graphs #################
###########################################################
##Figure 4
dframe10$High <- dframe10$high.protein
dframe10$High<- droplevels(dframe10$High)
levels(dframe10$High)<-c("Lower protein","Higher protein")
dframe10$High

#make data frame for annotation on one facet
dat_text <- data.frame(
  label = c("", "A"),
  High   = c("Higher protein", "Lower protein")
)

p1<-ggplot(dframe10, aes(x=Test)) +
  geom_bar(aes(fill = first.food.eaten ),position = "stack", width=0.75) +
  scale_fill_manual(values = c("#663300", "#E69F00", "#56B4E9", "#009E73")) +
  scale_x_discrete(name = "Test day",labels = c("Day 5","Day 10","Day 15","Day 35")) +
  facet_grid(~factor(High, levels= c("Higher protein", "Lower protein"))) +
  labs(x="Test day", y="Frequency") +
  theme_bw() + theme(text = element_text(size=22)) + theme(legend.position="none") +guides(fill="none") + 
  theme(strip.placement = "outside",strip.background =element_rect(fill=NA, colour = "white"),strip.text.x = element_text())+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black", linetype = "solid"))+
  geom_text(x = 4, y =20, aes(label = label), data = dat_text, size = 10)

p1

dframe110$High2 <- dframe110$secondhigh.protein
dframe110$High2
dframe110$High2<- droplevels(dframe110$High2)
levels(dframe110$High2)<-c("Lower protein","Higher protein")
dframe110$High2
#make data frame for annotation on one facet
dat_text2 <- data.frame(
  label = c("", "B"),
  High2   = c("Higher protein", "Lower protein")
)

p2<-ggplot(dframe110, aes(x=Test)) +
  geom_bar(aes(fill=second.food),position= "stack", width=0.75) +
  scale_fill_manual(values = c("#663300", "#E69F00", "#56B4E9", "#009E73")) +
  scale_x_discrete(name = "Test day",labels = c("Day 5","Day 10","Day 15","Day 35")) +
  scale_y_continuous(name=" ", limits=c(0, 20), breaks =seq(0,20,5)) +
  facet_grid(~factor(High2, levels= c("Higher protein", "Lower protein"))) +
  labs(x="Test day") +
  theme_bw() + theme(text = element_text(size=22)) + theme(legend.text=element_text(size=15),
    legend.title=element_text(size=15)) + guides(fill=guide_legend(title="Food choice"))  + 
  theme(strip.placement = "outside",strip.background =element_rect(fill=NA, colour = "white"),strip.text.x = element_text())+ 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black", linetype = "solid")) +
  geom_text(x = 4, y =20, aes(label = label), data = dat_text2, size = 10)
p2
#to graph in one plot
library(patchwork)
p1 + p2


###################### Statistical tests ########################################
#################################################################################

##### Do chicks prefer Marine over Terrestrial foods? #####
###########################################################
#Marine food =1
#Terrestrial =2 


###Exact binomial test for preference of Marine and Terrestrial food types
table(dframe10$Food.type)
binom.test(75,88, p= .5)

##Influence of diet treatment group, test day and individual identity on food preferences
###GLMM for food preference of Marine or Terrestrial food types

foodpref<- glmer(Food.type ~ Treatment + Test + 
                   (1|Chick.ID), data=dframe10,
                 family = binomial(link = "probit")) 
#singularity fit suggesting overfitting due to Chick ID
#Chick.ID controlling for repeat chick measures
#With Treatment group and Test trial 
#Not influential


summary(foodpref) 
AIC(foodpref) # logit= 85.025 cloglog= 85.102 probit =85.02 Probit offers best model fit
drop1(foodpref, test= "Chisq") #
par(mfrow=c(2,2))
plot(foodpref) 
check_model(foodpref) # 
exp(coef(foodpref)) #odds ratio back transformed
exp(confint(foodpref))


######## Individual foods preference within 4 food array offered ########
###Chi Squared goodness of fit
##All foods expected 25% likelihood of being eaten first
##Test days tested separately to avoid psuedo replicating chick.id
expected<- c(.25,.25,.25,.25)

#day 5
dframe6$Test<- droplevels(dframe6$Test)
five<- table(dframe6$first.food.eaten,dframe6$Test)
five
chisq.test(five, p=expected)
fisher.test(dframe6$first.food.eaten,dframe6$Test)

#day 10
ten<- table(dframe7$first.food.eaten)
ten
chisq.test(ten)

#day 15
fifteen<- table(dframe8$first.food.eaten)
fifteen
chisq.test(fifteen)

#day 35
thirtyfive<- table(dframe9$first.food.eaten)
thirtyfive
chisq.test(thirtyfive)

#This shows no difference in choice patterns across tests
dframe10$Test<- droplevels(dframe10$Test)
table0<- table(dframe10$Test,dframe10$first.food.eaten)
table0
chisq.test(table0)
fisher.test(table0)

### Post-hoc test to confirm which foods driving difference ###
#install.packages("RVAideMemoire")
library(RVAideMemoire)

chisq.theo.multcomp(five, p.method = "bonferroni")
chisq.theo.multcomp(ten, p.method = "bonferroni")
chisq.theo.multcomp(fifteen, p.method = "bonferroni")
chisq.theo.multcomp(thirtyfive, p.method = "bonferroni")


### Initial baseline test ###
dframe24$first.food.eaten<- droplevels(dframe24$first.food.eaten)
dframe24$first.food.eaten

expected<- c(.5,.5)
one<- table(dframe24$first.food.eaten)
one
chisq.test(one, p=expected)

#### Individual consistency in food preference ####
## Requires dataframe in 'long' form to compare individual chicks
## performance within the repeated food trials.
reshape <- read.csv("Rehab_Chicks Foodpreference_long.csv")
head(reshape) #puts an x in front of all numbers in column names
dim(reshape[1]) 
names(reshape)

# Remove Non Study birds
reshape <- reshape[reshape$Treatment!="Nonstud",]
dim(reshape[1]) #n= 27 study chicks

summary(reshape)
names(reshape)

#Changing variable classes to factors
reshape$Chick.ID <- as.factor(reshape$Chick.ID)
reshape$Treatment <- as.factor(reshape$Treatment)
reshape$X5first.food.eaten <- as.factor(reshape$first.food.eaten)
reshape$X5second.food <- as.factor(reshape$second.food)
reshape$X10first.food.eaten <- as.factor(reshape$X10first.food.eaten)
reshape$X10second.food <- as.factor(reshape$X10second.food)
reshape$X15first.food.eaten <- as.factor(reshape$X15first.food.eaten)
reshape$X15second.food <- as.factor(reshape$X15second.food)
reshape$X35first.food.eaten <- as.factor(reshape$X35first.food.eaten)
reshape$X35second.food <- as.factor(reshape$X35second.food)
reshape$highprotein5<- as.factor(reshape$highprotein5)
reshape$secondhigh.protein5 <- as.factor(reshape$secondhigh.protein5)
reshape$highprotein10<- as.factor(reshape$high.protein10)
reshape$secondhigh.protein10 <- as.factor(reshape$secondhigh.protein10)
reshape$highprotein15<- as.factor(reshape$high.protein15)
reshape$secondhigh.protein15 <- as.factor(reshape$secondhigh.protein15)
reshape$highprotein35<- as.factor(reshape$high.protein35)
reshape$secondhigh.protein35 <- as.factor(reshape$secondhigh.protein35)

#Mapping First and second food choice variables to recognise NAs and become numerics
levels(reshape$X5first.food.eaten)[levels(reshape$X5first.food.eaten)==""] <- NA
reshape$X5first.food.eaten
levels(reshape$X5second.food)[levels(reshape$X5second.food)==""] <- NA
reshape$X5second.food
levels(reshape$X10first.food.eaten)[levels(reshape$X10first.food.eaten)==""] <- NA
reshape$X10first.food.eaten
levels(reshape$X10second.food)[levels(reshape$X10second.food)==""] <- NA
reshape$X10second.food
levels(reshape$X15first.food.eaten)[levels(reshape$X15first.food.eaten)==""] <- NA
reshape$X15first.food.eaten
levels(reshape$X15second.food)[levels(reshape$X15second.food)==""] <- NA
reshape$X15second.food
levels(reshape$X35first.food.eaten)[levels(reshape$X35first.food.eaten)==""] <- NA
reshape$X35first.food.eaten
levels(reshape$X35second.food)[levels(reshape$X35second.food)==""] <- NA
reshape$X35second.food

mapping <- c("C" = 1, "F" = 2, "M" = 3, "B" = 4)
reshape$first5 <- mapping[reshape$X5first.food.eaten]
reshape$first10 <- mapping[reshape$X10first.food.eaten]
reshape$first15 <- mapping[reshape$X15first.food.eaten]
reshape$first35 <- mapping[reshape$X35first.food.eaten]
reshape$second5 <- mapping[reshape$X5second.food]
reshape$second10 <- mapping[reshape$X10second.food]
reshape$second15 <- mapping[reshape$X15second.food]
reshape$second35 <- mapping[reshape$X35second.food]

reshape$first5
reshape$X5first.food.eaten

### Interclass Correlation Coefficient ###
#Less than 0.50: Poor reliability
#Between 0.5 and 0.75: Moderate reliability
#Between 0.75 and 0.9: Good reliability
#Greater than 0.9: Excellent reliability

# First food choice
data1 <- data.frame(A=reshape$first5,
                   B=reshape$first10,
                   C=reshape$first15
                   )
data1

# Second food choice
data2 <- data.frame(A=reshape$second5,
                    B=reshape$second10,
                    C=reshape$second15
                    )
data2


#Comparing singular measurements of chick responses between trials
icc(data1, model="twoway", type = "consistency", unit = "single")
icc(data2, model="twoway", type = "consistency", unit = "single")

icc(data1, model="twoway", type = "agreement", unit = "single")
icc(data2, model="twoway", type = "agreement", unit = "single")


########################################################dframe10####################################################################################
########################## Do chicks prefer high protein to low protein options? ###
####################################################################################
#####Binomial tests comparing test days vs preference for high protein y/n

##First choice
first<-dframe10$high.protein
first


#day5
binom.test(20,25, p=.5)
#day10
binom.test(17,26, p=.5)
#day15
binom.test(13,24, p=.5)
#day35
binom.test(10,13, p=.5)

#Second choice
second<-dframe10$secondhigh.protein
second
table (first,second, dframe10$Test)


#second and first choice combined
#day5
binom.test(25,25, p=.83)
#day10
binom.test(26,26, p=.83)
#day15
binom.test(23,24, p=.83)
#day35
binom.test(12,13, p=.83)


########################################################################
###### Interobserver Comparison ######

dframe1 <- read.csv("Observer comparison.csv", header = T)
summary(dframe1)
head(dframe1)
names(dframe1)

#### Remove Non Study gulls from list (no data there)
dframe1 <- dframe1[dframe1$Treatment!="Nonstud",]
dim(dframe1[1])

###Subset data to trials without Nas for food choice
dframe2 <- dframe1[!is.na(dframe1$first.food.eaten),]
summary(dframe2)
dim(dframe2[1])


#change the class of variables
dframe1$emma.food.latency <- as.integer(dframe1$emma.food.latency) 
dframe1$frankie.food.latency <- as.integer(dframe1$frankie.food.latency) 

#making the first food eaten variable numeric
mapping <- c("F" = 1, "M" = 2, "C"= 3, "B"= 4)
dframe1$emma.food.numeric <- mapping[dframe1$first.food.eaten]
dframe1$frankie.food.numeric <- mapping[dframe1$frankie.first.choice]

####ggplots####
library(ggplot2)

length(dframe2)
ggplot(dframe2) +
  geom_point(aes(x = emma.food.latency, y=frankie.food.latency)) + 
  theme(axis.text.x=element_text(angle=0,hjust=1)) +
  geom_abline()+labs(title="Scores of gull latency to approach food", x="Emma", y="Frankie")

#Food item choice
ggplot(dframe1) +
  geom_point(aes(x = emma.food.numeric, y=frankie.food.numeric)) + 
  theme(axis.text.x=element_text(angle=0,hjust=1)) +
  geom_abline() + labs(title="First food eaten", x="Emma", y="Frankie")


################ Comparison between raters ###################
library("irr")


###Interclass Correlation Coefficient###
#Less than 0.50: Poor reliability
#Between 0.5 and 0.75: Moderate reliability
#Between 0.75 and 0.9: Good reliability
#Greater than 0.9: Excellent reliability


#Frankie and Emma
data <- data.frame(A=dframe2$emma.food.latency,
                   B=dframe2$frankie.food.latency)
data

data1 <- data.frame(A=dframe1$emma.food.numeric,
                    B=dframe1$frankie.food.numeric)
data1


#Comparing singular measurements of chick responses between trials
icc(data1, model="twoway", type = "agreement", unit = "single")
# Excellent reliability and absolute agreement between raters are observed using two-way 
# random effect models for each measure with ICC coefficient of 0.917 (95% Cl 0.884-0.941) p = 1.01e-52 
# First food choice excellent reliability between raters ICC 0.921 (95% Cl 0.887 - 0.945)  p = 5.06e-47

################################################################################
#################### Chick Biometric comparisons ###############################

#### Chick weight ####
dframe1 <- read.csv("REHAB_Chicksweight.csv")
summary(dframe1)
head(dframe1)
names(dframe1)
#%>% #make by control shift M

#change the class of variables
dframe1[,3:81] <- lapply(dframe1[,3:81], as.integer)

#Remove nonstudy chicks and empty rows
dframe1<-dframe1[dframe1$chick.ID!="" & dframe1$chick.ID!="NONSTUD1" & dframe1$chick.ID!="NONSTUD2",]
dframe1$chick.ID

#####Testing weight differences between groups ######
marine<- dframe1[dframe1$Group=="Marine",]
terra<- dframe1[dframe1$Group=="Terra",]

hist(marine$last.weight)
hist(terra$last.weight)

marine$last.weight<-as.numeric(marine$last.weight)
terra$last.weight<-as.numeric(terra$last.weight)

#wilcoxon test n = all 27 study chicks 14 marine, 13 terra
wilcox.test(marine$last.weight,terra$last.weight, exact=TRUE, conf.int = TRUE)

#test for chicks 15 days old n = 15, 8 marine and 7 terra
wilcox.test(marine$X15,terra$X15, exact=TRUE, conf.int = TRUE)

#Test for chicks 5 days old n = 8, 4 from both study groups
wilcox.test(marine$X5,terra$X5, exact=TRUE, conf.int = TRUE)

###### To make weight graph ######
####Make dataset long from wide
library("reshape2")
dframe1b<-dframe1[dframe1$chick.ID!="",]
dframe2<-melt(dframe1b, id.vars = "chick.ID", measure.vars = c("X1","X2","X3",
                                                               "X4","X5","X6","X7",
                                                               "X8","X9","X10","X11",
                                                               "X12","X13","X14","X15",
                                                               "X16","X17","X18","X19",
                                                               "X20","X21","X22","X23",
                                                               "X24","X25","X26","X27",
                                                               "X28","X29","X30","X31",
                                                               "X32","X33","X34","X35",
                                                               "X36","X37","X38","X39",
                                                               "X40","X41","X42","X43",
                                                               "X44","X45","X46","X47",
                                                               "X48","X49","X50","X51",
                                                               "X52","X53","X54","X55",
                                                               "X56","X57","X58","X59",
                                                               "X60","X61","X62","X63",
                                                               "X64","X65","X66","X67",
                                                               "X68","X69","X70","X71",
                                                               "X72","X73","X74","X75",
                                                               "X76","X77","X78","X79"))

#Age<-dframe2$variable
#Weight<-dframe2$value
names(dframe2)[names(dframe2)=="variable"] <- "Age"
names(dframe2)[names(dframe2)=="value"] <- "Weight"

#Ordering data by chick and age
dframe2 <- dframe2[order(dframe2$chick.ID, dframe2$Age), ]
head(dframe2)
summary(dframe2)

#Changing name of levels in age
# Across all columns, replace all instances of "X" with " "
levels(dframe2$Age) <- gsub("X", " ", levels(dframe2$Age))
dframe2$Age
class(dframe2$Age) #still factor

#Making age integer
dframe2$Age<-as.integer(dframe2$Age)

###Graphing###
library(ggplot2)
# Basic line graph
ggplot(data=dframe2, aes(x=Age, y=Weight, group=chick.ID, color=chick.ID)) +
  geom_line()
## This would have the same result as above
# ggplot(data=dframe2, aes(x=Weight, y=Age)) +
#     geom_line(aes(group=chick.ID))

# Add points
ggplot(data=dframe2, aes(x=Age, y=Weight, group=chick.ID, color=chick.ID)) +
  geom_line() +
  geom_point(pch="+", cex=4) +
  xlab("Chick age (days)") + ylab("Chick weight (g)") +
  scale_y_continuous(limits=c(0,1250)) + 
  scale_x_continuous(breaks=seq(0,85,10),limits=c(0,85))  +
  theme_bw() + theme(strip.background =element_rect(fill="gray100")) + 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(legend.position="none")     # This must go after theme_bw


## separate the chicks by group
#merge two spreadsheets
food<-  read.csv("REHAB_Chicks Foodpreference_reshape.csv")
names(food)

names(dframe2)
dframe2$chick.ID
dframe2$Treatment<-NA
for(i in 1:dim(dframe2)[1]) {
  
  focal.chick<-dframe2$chick.ID[i]
  
  temp.chick<-food[food$Chick.ID==focal.chick,]
  
  if(dim(temp.chick)[1]>0) {
    
    dframe2$Treatment[i]<-temp.chick$Treatment[1]
  }
  
  print(i)
}

summary(dframe2) 
temp.chick$Treatment
focal.chick

ggplot(data=dframe2, aes(x=Age, y=Weight, group=chick.ID, color=chick.ID)) +
  geom_line() +
  geom_point(pch="+", cex=4) +
  xlab("Chick age (days)") + ylab("Chick weight (g)") +
  scale_y_continuous(limits=c(0,1250)) + 
  scale_x_continuous(breaks=seq(0,85,10),limits=c(0,85))  +
  theme_bw() + theme(strip.background =element_rect(fill="gray100")) + 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(legend.position="none") 



##Colour blind palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
#To use for fills, add
scale_fill_manual(values=cbPalette)
# To use for line and point colors, add
scale_colour_manual(values=cbPalette)

#########################################################################################################
#### Tarsus ####
dframe3 <- read.csv("REHAB_Chickstarsus.csv")
summary(dframe3)
head(dframe3)
names(dframe3)
#names(dframe3)[names(dframe3)=="X"] <- "chick.ID"

#change the class of variables
#dframe3[,2:80] <- lapply(dframe3[,2:80], as.numeric)

#Remove nonstudy chicks and empty rows
dframe3<-dframe3[dframe3$Chick.ID!="" & dframe3$Chick.ID!="NONSTUD1" & dframe3$Chick.ID!="NONSTUD2",]
dframe3$Chick.ID

#####Testing weight differences between groups ######
dframe3$Group


marine<- dframe3[dframe3$Group=="Marine",]
terra<- dframe3[dframe3$Group=="Terra",]

hist(marine$last.tar)
hist(terra$last.tar)


#wilcoxon test last chick tarsus measurement
#chick age range
summary(dframe3$last.age)

wilcox.test(marine$last.tar,terra$last.tar, exact=TRUE, conf.int = TRUE)

#chick ages 14-16 days comparison
wilcox.test(marine$tar.age14.16,terra$tar.age14.16, exact=TRUE, conf.int = TRUE)



###### TO make Tarsus graph ######
####Make dataset long from wide
dframe4<-melt(dframe3, id.vars = "chick.ID", measure.vars = c("X1","X2","X3","X4"
                                                              ,"X5","X6","X7","X8",
                                                              "X9","X10","X11","X12",
                                                              "X13","X14","X15","X16",
                                                              "X17","X18","X19","X20",
                                                              "X21","X22","X23","X24",
                                                              "X25","X26","X27","X28",
                                                              "X29","X30","X31","X32",
                                                              "X33","X34","X35","X36",
                                                              "X37","X38","X39","X40",
                                                              "X41","X42","X43","X44",
                                                              "X45","X46","X47","X48",
                                                              "X49","X50","X51","X52",
                                                              "X53","X54","X55","X56",
                                                              "X57","X58","X59","X60",
                                                              "X61","X62","X63","X64",
                                                              "X65","X66","X67","X68",
                                                              "X69","X70","X71","X72",
                                                              "X73","X74","X75","X76",
                                                              "X77","X78","X79"))

names(dframe4)[names(dframe4)=="variable"] <- "Age"
names(dframe4)[names(dframe4)=="value"] <- "Tarsus"

#Ordering data by chick and age
dframe4 <- dframe4[order(dframe4$chick.ID, dframe4$Age), ]
head(dframe4)

#Changing name of levels in age
# Across all columns, replace all instances of "X" with " "
levels(dframe4$Age) <- gsub("X", " ", levels(dframe4$Age))
dframe4$Age
class(dframe4$Age) #still factor

#Making age integer
dframe4$Age<-as.integer(dframe4$Age)


###Graph###
ggplot(data=dframe4, aes(x=Age, y=Tarsus, group=chick.ID, color=chick.ID)) +
  geom_line() + 
  geom_point(pch="+", cex=3) +
  xlab("Chick age (days)") + ylab("Minimal tarsus measurement (mm)") +
  scale_y_continuous(limits=c(0,75)) + 
  scale_x_continuous(breaks=seq(0,60,10),limits=c(0,60))  +
  theme_bw() + theme(strip.background =element_rect(fill="gray100")) + 
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  theme(legend.position="none")     # This must go after theme_bw





#########################################################################################################
########################## END ####################################