install.packages("car")
library(robustHD)
library(robustbase)
library(MASS)
library(car)
library(ggplot2)
install.packages("emmeans")
library(emmeans)
library(multcomp)
library(Rmisc)
library("ggpubr")
#######################################################################
#### adult nutrition experiment ####
#######################################################################

#####################################################
####   weight analysis and plots   ####
#####################################################

#attaching the file
hangry.data<-read.csv(file.choose())
names(hangry.data)

#turning food deprivation into a factor
hangry.data$food.deprivation.f<-as.factor(hangry.data$food.deprivation)

#### wet weight ####
hangry.data$av.wet.weight<-(hangry.data$Fly.1.wet.mass+hangry.data$Fly.2.wet.mass)/2

#looking at the spread of the data
hist(hangry.data$av.wet.weight)
boxplot(hangry.data$av.wet.weight, main="wet weight")

#wet weight against food deprivation 
wetweight.model.factor<-lm(av.wet.weight~food.deprivation.f+block, data=hangry.data)
summary(wetweight.model.factor)
Anova(wetweight.model.factor, type="II")
#checking distribution of residuals
hist(wetweight.model.factor$residuals)
qqnorm(residuals(wetweight.model.factor))
qqline(residuals(wetweight.model.factor))

emmeans(wetweight.model.factor, list(pairwise ~ food.deprivation.f), adjust = "tukey")

#dry weight
hangry.data$av.dry.weight<-(hangry.data$Fly.1.dry.mass+hangry.data$Fly.2.dry.mass)/2
hist(hangry.data$av.dry.weight)
boxplot(hangry.data$av.dry.weight, main="dry weight")

#dry weight against food deprivation
dryweight.model.factor<-lm(av.dry.weight~food.deprivation.f+block, data=hangry.data)
summary(dryweight.model.factor)
Anova(dryweight.model.factor, type="II")
#checking distribution of residuals
hist(dryweight.model.factor$residuals)
qqnorm(residuals(dryweight.model.factor))
qqline(residuals(dryweight.model.factor))

emmeans(dryweight.model.factor, list(pairwise ~ food.deprivation.f), adjust = "tukey")

########################################################################
#### behavior ####
########################################################################

#making a variable that includes all forms of aggression
hangry.data$aggression<-hangry.data$lunge+hangry.data$chase+hangry.data$fence+hangry.data$tussle

#turning it into a rate per minute
hangry.data$aggression.rate<-(hangry.data$aggression/(hangry.data$number.of.scans*3)*60)

#turning time on food patch into rate per minute
hangry.data$food.rate<-(hangry.data$feeding/(hangry.data$number.of.scans*3)*60)

#####################################################
#### aggression rate plots and analysis ####
#####################################################

#looking at the spread of the data
hist(hangry.data$aggression.rate)
adjbox(hangry.data$aggression.rate)

#aggression against food deprivation
hangry.data$food.deprivation.f<-as.factor(hangry.data$food.deprivation)
full.model.factorlm<-lm(sqrt(aggression.rate)~food.deprivation.f+block,  data=hangry.data)

summary(full.model.factorlm)
Anova(full.model.factorlm)
#food deprivation and block are signficant
#checking distribution of residuals
par(mfrow=c(2,2))
plot(full.model.factorlm)
par(mfrow=c(1,1))
hist(residuals(full.model.factorlm))
qqnorm(residuals(full.model.factorlm))
qqline(residuals(full.model.factorlm))
summary(full.model.factorlm$residuals)

emmeans(full.model.factorlm, list(pairwise ~ food.deprivation.f), adjust = "tukey")


#### is this due to weight? ####

wet.weight.model<-lm(sqrt(aggression.rate)~av.wet.weight+block,  data=hangry.data)
dry.weight.model<-lm(sqrt(aggression.rate)~av.dry.weight+block,  data=hangry.data)

#does wet weight predict aggression? 
summary(wet.weight.model)
Anova(wet.weight.model)
hist(residuals(wet.weight.model))
qqnorm(residuals(wet.weight.model))
qqline(residuals(wet.weight.model))

#does dry weight predict aggression? n
summary(dry.weight.model)
Anova(dry.weight.model)
hist(residuals(dry.weight.model))
qqnorm(residuals(dry.weight.model))
qqline(residuals(dry.weight.model))

#are there effects of dry weight over and above the nutrient manipulation?
full.model.withdw<-lm(sqrt(aggression.rate)~food.deprivation.f+av.dry.weight+block,  data=hangry.data)
summary(full.model.withdw)
anova(full.model.withdw)
#yes

library(jtools)
effect_plot(full.model.withdw, pred =av.dry.weight, interval = TRUE, plot.points = TRUE)


#### lunges ####
#did they lunge or not? 
hangry.data$LungeOrNot<-ifelse(hangry.data$lunge==0,0,1)
#the variable 'LungeOrNot' is a 0 or 1 to say if they lunged or not
# this is binomial 
summary(as.factor(hangry.data$LungeOrNot))

LungeOrNot2 <- glm(LungeOrNot ~food.deprivation.f2+block,family=binomial(link='logit'),data=hangry.data)
summary(LungeOrNot2)
Anova(LungeOrNot2)
emmeans(LungeOrNot2, list(pairwise ~ food.deprivation.f), adjust = "tukey")


#fencing
#did they fence or not? 
#making a fence or not score
hangry.data$FenceOrNot<-ifelse(hangry.data$fence==0,0,1)
# this is binomial 
summary(as.factor(hangry.data$FenceOrNot))

FenceOrNot2 <- glm(FenceOrNot ~food.deprivation.f2+block,family=binomial(link='logit'),data=hangry.data)
summary(FenceOrNot2)
Anova(FenceOrNot2)

emmeans(FenceOrNot2, list(pairwise ~ food.deprivation.f), adjust = "tukey")

#chasing 
#did they chase or not? 
#making a chase or not score
hangry.data$ChaseOrNot<-ifelse(hangry.data$chase==0,0,1)
#the variable 'LungeOrNot' is a 0 or 1 to say if they lunged or not
# this is binomial 
summary(as.factor(hangry.data$ChaseOrNot))

ChaseOrNot2 <- glm(ChaseOrNot ~food.deprivation.f2+block,family=binomial(link='logit'),data=hangry.data)
summary(ChaseOrNot2)
Anova(ChaseOrNot2)

emmeans(ChaseOrNot2, list(pairwise ~ food.deprivation.f), adjust = "tukey")


hangry.data$tussle.rate<-(hangry.data$tussle/(hangry.data$number.of.scans*3)*60)
hangry.data$TussleOrNot<-ifelse(hangry.data$tussle.rate==0,0,1)
# this is binomial 
summary(as.factor(hangry.data$TussleOrNot))
#only pairs 2 tussled

#were each of the aggressive behaviours correlated?
chisq.test(hangry.data$LungeOrNot, hangry.data$FenceOrNot)
chisq.test(hangry.data$LungeOrNot, hangry.data$ChaseOrNot)
chisq.test(hangry.data$ChaseOrNot, hangry.data$FenceOrNot)

#####################################################
#### plotting food patch against food deprivation ####
#####################################################

#looking at the spread of the data
hist(hangry.data$food.rate)
boxplot(hangry.data$food.rate)
adjbox(hangry.data$food.rate)

fullfood.model.factor<-lm(sqrt(food.rate)~food.deprivation.f+block, data=hangry.data)
summary(fullfood.model.factor)
Anova(fullfood.model.factor)
#food deprivation and block are signficant
#checking distribution of residuals
hist(residuals(fullfood.model.factor))
qqnorm(residuals(fullfood.model.factor))
qqline(residuals(fullfood.model.factor))
emmeans(fullfood.model.factor, "food.deprivation.f")
emmeans(fullfood.model.factor, list(pairwise ~ food.deprivation.f), adjust = "tukey")


#### are there effects of time on food patch above and beyond nutrition? ####
fullfood.model.factor<-lm(sqrt(aggression.rate)~food.deprivation.f+food.rate+block, data=hangry.data)
summary(fullfood.model.factor)
anova(fullfood.model.factor)


#correlation between time on food and aggression

cor.test(hangry.data$aggression.rate,hangry.data$food.rate, method="spearman")
shapiro.test(hangry.data$food.rate)
shapiro.test(hangry.data$aggression.rate)


#### plots ####

hangry.data$food.deprivation.f2<-hangry.data$food.deprivation.f
hangry.data$food.deprivation.f2 <- factor(hangry.data$food.deprivation.f2, levels = c("0", "24", "48","72","120"),
                                          labels = c("0", "24", "48","72","120+"))


Aggression.rate.plot<-ggplot(hangry.data, aes(x = food.deprivation.f2, y = aggression.rate))+labs(y="aggression.rate.adj", x= "food deprivation")+theme_classic()+
  geom_jitter(shape=16,col="grey60", position=position_jitter(0.2))+
  labs(y="Aggression Rate", x="Food Deprivation Duration (h)")+
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 35, size = 7, colour = "grey30")+
  annotate("text", label = "ac", x = 2, y = 35, size = 7, colour = "grey30")+
  annotate("text", label = "bc", x = 3, y = 35, size = 7, colour = "grey30")+
  annotate("text", label = "bd", x = 4, y = 35, size = 7, colour = "grey30")+
  annotate("text", label = "bc", x = 5, y = 35, size = 7, colour = "grey30")

lunge.prob.plot<-ggplot(LungeOrNot2, aes(x = food.deprivation.f2, y = LungeOrNot))+theme_classic()+
  labs(y="Probability of Lunging", x="Food Deprivation Duration (h)")+
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  scale_y_continuous(
    labels = scales::number_format(accuracy = 0.1))+
  annotate("text", label = "a", x = 1, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "ab", x = 2, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "ab", x = 3, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 4, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "ab", x = 5, y = 1, size = 7, colour = "grey30")

fence.prob.plot<-ggplot(FenceOrNot2, aes(x = food.deprivation.f2, y = FenceOrNot))+theme_classic()+
  labs(y="Probability of Fencing", x="Food Deprivation Duration (h)")+
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  scale_y_continuous(
    labels = scales::number_format(accuracy = 0.1))+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "ab", x = 2, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 3, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 4, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "ab", x = 5, y = 1, size = 7, colour = "grey30")


chase.prob.plot<-ggplot(ChaseOrNot2, aes(x = food.deprivation.f2, y = ChaseOrNot))+theme_classic()+
  labs(y="Probability of Chasing", x="Food Deprivation Duration (h)")+
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  scale_y_continuous(
    labels = scales::number_format(accuracy = 0.1))+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "a", x = 2, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "a", x = 3, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "a", x = 4, y = 1, size = 7, colour = "grey30")+
  annotate("text", label = "a", x = 5, y = 1, size = 7, colour = "grey30")

food.adult.plot<-ggplot(hangry.data, aes(x = food.deprivation.f2, y = food.rate))+theme_classic()+
  geom_jitter(shape=16,col="grey60", position=position_jitter(0.2))+
  labs(y="Food Patch Occupancy", x="Food Deprivation Duration (h)")+
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 26, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 2, y = 26, size = 7, colour = "grey30")+
  annotate("text", label = "bd", x = 3, y = 26, size = 7, colour = "grey30")+
  annotate("text", label = "cd", x = 4, y = 26, size = 7, colour = "grey30")+
  annotate("text", label = "d", x = 5, y = 26, size = 7, colour = "grey30")

food.cor.plot<-ggplot(hangry.data, aes(x=food.rate, y=aggression.rate))   +labs(y="Aggression Rate", x="Food Patch Occupancy") +  theme_classic()+
  geom_point(col="grey50")+stat_smooth(method = "lm", col = "red")+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )

ggarrange(Aggression.rate.plot,lunge.prob.plot,fence.prob.plot,chase.prob.plot, labels=c("A","B","C","D"), nrow=2,ncol=2)

ggarrange(Aggression.rate.plot,lunge.prob.plot,fence.prob.plot,chase.prob.plot,food.adult.plot,food.cor.plot, labels=c("A","B","C","D","E","F"), nrow=3,ncol=2)


wetweightplot<-ggplot(wetweight.model.factor, aes(x=food.deprivation.f, y=av.wet.weight))  +labs(y="Wet Mass (mg)", x="Food Deprivation Duration (hrs)")+theme_classic()+
  geom_violin()+
  geom_jitter(shape=16, position=position_jitter(0.2), col="grey50")+
  scale_x_discrete(name ="Food Deprivation Duration",
                   labels=c("0","24","48","72","120+")) +
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 1.2, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 2, y = 1.2, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 3, y = 1.2, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 4, y = 1.2, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 5, y = 1.2, size = 7, colour = "grey30")

dryweightplot<-ggplot(dryweight.model.factor, aes(x=food.deprivation.f, y=av.dry.weight))  +labs(y="Dry Mass (mg)", x="Food Deprivation Duration (hrs)")+theme_classic()+
  geom_violin()+
  geom_jitter(shape=16, position=position_jitter(0.2),col="grey50")+
  scale_color_grey()+
  scale_x_discrete(name ="Food Deprivation Duration",
                   labels=c("0","24","48","72","120+")) +
  stat_summary(fun.y=mean, geom="point", shape=20, size=5, col="black")+
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", col="black", size=0.9,width=0.1)+
  theme(axis.text=element_text(size=rel(1.4)),axis.title=element_text(size=rel(1.4)) )+
  annotate("text", label = "a", x = 1, y = 0.40, size = 7, colour = "grey30")+
  annotate("text", label = "b", x = 2, y = 0.40, size = 7, colour = "grey30")+
  annotate("text", label = "c", x = 3, y = 0.40, size = 7, colour = "grey30")+
  annotate("text", label = "d", x = 4, y = 0.40, size = 7, colour = "grey30")+
  annotate("text", label = "d", x = 5, y = 0.40, size = 7, colour = "grey30")

ggarrange(wetweightplot,dryweightplot, labels=c("A","B"), nrow=1,ncol=2)


cor.56<-ggplot(hangry.data, aes(x=av.wet.weight, y=aggression.rate))  +  theme_classic()+labs(x="Average Pair Wet Mass (mg)", y="Aggression Rate")+
  geom_point(col="grey50")+
  theme(axis.text=element_text(size=rel(1.2)),axis.title=element_text(size=rel(1.4)) )
cor.95<-ggplot(hangry.data, aes(x=av.dry.weight, y=aggression.rate))  +  theme_classic()+labs(x="Average Pair Dry Mass (mg)", y="Aggression Rate")+
  geom_point(col="grey50")+stat_smooth(method = "lm", formula = y ~ x, col = "black")+
  theme(axis.text=element_text(size=rel(1.2)),axis.title=element_text(size=rel(1.4)) )

ggarrange(cor.56,cor.95, labels=c("A","B"), nrow=1,ncol=2)