#### R CODE FOR IOANNOU ET AL. COPYRIGHT 2025, CHRISTOS IOANNOU, ALL RIGHTS RESERVED.

#### LIBRARIES

library(MASS) 
library(stats) 
library(cluster) 
library(graphics)
# install.packages("lme4" ,lib="FILE PATH/R packages")
library(lme4,lib.loc="FILE PATH/R packages")
# install.packages("glmmTMB",lib="FILE PATH/R packages",dependencies=T)
library(glmmTMB,lib.loc="FILE PATH/R packages")
# install.packages("gap.datasets",lib="FILE PATH/R packages",dependencies=T)
library(gap.datasets,lib.loc="FILE PATH/R packages")
# install.packages("gap",lib="FILE PATH/R packages",dependencies=T)
library(gap,lib.loc="FILE PATH/R packages")
# install.packages("DHARMa",lib="FILE PATH/R packages",dependencies=T)
library(DHARMa,lib.loc="FILE PATH/R packages")
# install.packages("bbmle",lib="FILE PATH/R packages",dependencies=T)
library(bbmle,lib.loc="FILE PATH/R packages")
# install.packages("lmerTest",lib="FILE PATH/R packages",dependencies=T)
library(lmerTest,lib.loc="FILE PATH/R packages")
# install.packages("boot",lib="FILE PATH/R packages",dependencies=T)
library(boot,lib.loc="FILE PATH/R packages")


#### STATISTICAL ANALYSIS

Fdata <- read.csv(file="SupData.csv",stringsAsFactors=T)
nrow(Fdata[!complete.cases(Fdata),])
	# 19 CASES INCOMPLETE
str(Fdata)
length(unique(Fdata$Trial_number))
table(Fdata$Trial_number)
	# 120 TRIALS, ALTHOUGTH ONLY ONE PRESENTATION IN ONE TRIAL

Fdata$Min_dist_stimulus_cm <- Fdata$Min_dist_stimulus * 136 / 1690.157
Fdata$Med_distance_group_centroid_cm <- Fdata$Med_distance_group_centroid * 136 / 1690.157
	# CONVERSION TO REAL-WORLD DIMS: 1690.157 PIXELS FOR LENGTH OF ARENA, 136 CM IS REAL-WORLD LENGTH
Fdata$Med_speed_of_first_to_reach_cm <- Fdata$Med_speed_of_first_to_reach * 136 / 1690.157
	# SPEED IS CM PER SECOND
plot(Fdata$Min_dist_stimulus_cm~Fdata$Min_dist_stimulus)
plot(Fdata$Med_distance_group_centroid_cm~Fdata$Med_distance_group_centroid)
	# CHECK
hist(Fdata$Med_speed_of_first_to_reach_cm)


## DO THE EXPERIMENTAL VARIABLES (MANIPULATED OR MEASURED) AFFECT LATENCY TO ATTACK?

length(unique(Fdata$Trial_number))
table(Fdata$Trial_number)

m1<-glmer.nb(Latency_of_attack~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
simulationOutput <- simulateResiduals(m1,n=1000)
testDispersion(simulationOutput)
	# OVERDISPERSED WITH NEG BIN

m1<-lmer(log(Latency_of_attack)~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
simulationOutput <- simulateResiduals(m1,n=1000)
plot(simulationOutput)
	# ASSUMPTIONS MET
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
m3<-lmer(log(Latency_of_attack)~Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
m4<-lmer(log(Latency_of_attack)~          scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
m5<-lmer(log(Latency_of_attack)~                            scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
ICtab(m1,m2,m3,m4,m5,type=c("AICc"))
	# BOTH MAIN EFFECTS IMPORTANT, MODEL WITH INTERACTION IS CLOSE, BUT DOESN'T IMPROVE MODEL VS THE MAIN EFFECTS ONLY
   dAICc df
m2  0.0  9 
m1  0.9  10
m4 32.5  8 
m3 67.4  8 
m5 86.7  7 
length(residuals(m2))

summary(m2)
confint(m2)


## IS THE EFFECT OF GROUP SIZE EXPLAINED BY BEARING AND DISTANCE JUST BEFORE THE PRESENTATION?

m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# JUST GROUP SIZE
m2a<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)+Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH BEHAVIOURAL VARIABLE AS ANOTHER MAIN EFFECT
m2b<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)+Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH Min_dist_stimulus AND WITHOUT GROUP SIZE
m2c<-lmer(log(Latency_of_attack)~scale(Min_bearing_to_stimulus)+Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH Min_bearing_to_stimulus AS ANOTHER MAIN EFFECT
m2d<-lmer(log(Latency_of_attack)~scale(Min_bearing_to_stimulus)+Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH Min_bearing_to_stimulus AND WITHOUT GROUP SIZE
m2e<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)+scale(Min_bearing_to_stimulus)+Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH Min_bearing_to_stimulus AND Min_dist_stimulus AS MAIN EFFECTS
simulationOutput <- simulateResiduals(m2e,n=1000)
plot(simulationOutput)
m2f<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)+scale(Min_bearing_to_stimulus)+Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# WITH Min_bearing_to_stimulus AND Min_dist_stimulus AS MAIN EFFECTS AND WITHOUT GROUP SIZE
ICtab(m2,m2a,m2b,m2c,m2d,m2e,m2f,type=c("AICc"))
	# ADDING Min_dist_stimulus AND Min_dist_stimulus IMPROVES THE MODEL, BUT THEN REMOVING GROUP SIZE STILL HAS A HUGE NEGATIVE EFFECT
    dAICc df
m2e  0.0  11
m2c 10.1  10
m2a 39.7  10
m2f 41.3  10
m2  53.3  9 
m2d 63.6  9 
m2b 92.5  9 
length(residuals(m2e))

summary(m2e)
confint(m2e)



## DO THE BEHAVIOURAL VARIABLES INTERACT WITH TURBIDITY AND GROUP SIZE?

# Med_distance_group_centroid
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Med_distance_group_centroid_cm)*Treatment+scale(Med_distance_group_centroid_cm)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Med_distance_group_centroid_cm)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Med_distance_group_centroid_cm)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Med_distance_group_centroid_cm)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Med_distance_group_centroid_cm)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3 WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
# plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# BOTH INTERACTIONS IMPROVE THE MODEL FIT
    dAICc df
m21  0.0  12
m25  4.1  14
m23  4.1  11
m22 10.2  11
m2  11.6  9 
m24 13.7  10
length(residuals(m2))

summary(m21)
confint(m21)


# FIG. 2

tiff("Fig 2.tiff", width=8.7,height=11, units = "cm", type="windows", compression="lzw", res=600)
par(oma=c(0, 3, 0, 0))
par(mar=c(3, 1, 1.3, 0.4))
layout( matrix(c(1,1,1,2,3,4,5,6,7),nrow=3,byrow=T) )

boxplot(Fdata$Latency_of_attack~Fdata$Treatment*Fdata$Group_size,las=1,log="y",xlab="",ylab="",xaxt="n",yaxt="n"
	,col=rep(c(rgb(0/255,204/255,204/255,alpha=0),rgb(128/255,128/255,128/255,alpha=0.4)),3)
	,border=rep(c(rgb(0/255,204/255,204/255),"black"),3))
axis(side=1,at=1:6,labels=rep(c("Clear","Turbid"),3),mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=1,at=c(1.5,3.5,5.5),labels=c("2 fish","6 fish","10 fish"),mgp=c(1.5,0.6,0),las=1,tcl=0,line=1,lty="blank")
axis(side=2,labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("A",side=3,line=0.07,cex=1,adj=0,font=2)

fixef(m21)
Behaviour <- seq(from = min(scale(Fdata$Med_distance_group_centroid_cm),na.rm=T), to = max(scale(Fdata$Med_distance_group_centroid_cm),na.rm=T), length=100)
	# STILL WITH scale()
Turbid <- c(0,1)
Scaled_group_sizes <- sort(unique(scale(Fdata$Group_size)))
GS <- Scaled_group_sizes

Plot_data <- expand.grid(Behaviour,Turbid,GS)
colnames(Plot_data) <- c("Behaviour","Turbid","GS")

Predicted <- fixef(m21)[9]*Plot_data$Behaviour*Plot_data$GS +
	fixef(m21)[8]*Plot_data$Behaviour*Plot_data$Turbid +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*Plot_data$GS +
	fixef(m21)[3]*Plot_data$Turbid +
	fixef(m21)[2]*Plot_data$Behaviour +
	fixef(m21)[1]	
range(predict(m21))
range(Predicted)
	# CHECK: CLOSE MATCH
range(log(Fdata$Latency_of_attack))

Plot_data$Behaviour_not_scaled <- rep( seq(from = min(Fdata$Med_distance_group_centroid_cm,na.rm=T), to = max(Fdata$Med_distance_group_centroid_cm,na.rm=T), length=100) ,6)
	# FOR PLOTTING ONLY
# plot(Plot_data$Behaviour_not_scaled~Plot_data$Behaviour)
	# CHECK: IS A LINEAR SCALING

mtext("Latency to attack (s)",side=2,line=1.5,adj=0.5,outer=T,cex=0.8)

par(mar=c(3, 1, 1.3, 0.4))
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=range(Fdata$Med_distance_group_centroid_cm,na.rm=T)
	,xlab="",ylab="",main="2 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==2]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Clear" & Fdata$Group_size==2]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[1]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[1]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[1]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[1]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,20,40,60),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("B",side=3,line=0.07,cex=1,adj=0,font=2)

par(mar=c(3, 1, 1.3, 0.4))		
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=range(Fdata$Med_distance_group_centroid_cm,na.rm=T)
	,xlab="",ylab="",main="6 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==6]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Clear" & Fdata$Group_size==6]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[2]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[2]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[2]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[2]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,20,40,60),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
par("cex")
mtext("Median distance to group centroid (cm)",side=1,line=1.7,adj=0.5,outer=F,cex=0.8)

par(mar=c(3, 1, 1.3, 0.4))
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==10]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Turbid" & Fdata$Group_size==10]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=range(Fdata$Med_distance_group_centroid_cm,na.rm=T)
	,xlab="",ylab="",main="10 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==10]~Fdata$Med_distance_group_centroid_cm[Fdata$Treatment=="Clear" & Fdata$Group_size==10]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[3]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[3]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[3]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[3]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,20,40,60),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
	# CONTINUE PLOT AFTER POLARISATION ANALYSIS NEXT


# Med_group_polarisation
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Med_group_polarisation)*Treatment+scale(Med_group_polarisation)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Med_group_polarisation)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Med_group_polarisation)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Med_group_polarisation)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Med_group_polarisation)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3-WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
# plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# BOTH INTERACTIONS IMPROVE THE MODEL FIT
	# THREE-WAY INTERACTION DOESN'T IMPROVE FIT OVER THIS
     dAICc df
m21  0.0  12
m25  3.8  14
m22  9.9  11
m23 10.9  11
m24 20.5  10
m2  23.3  9 
length(residuals(m2))

summary(m21)
confint(m21)

fixef(m21)
Behaviour <- seq(from = min(scale(Fdata$Med_group_polarisation),na.rm=T), to = max(scale(Fdata$Med_group_polarisation),na.rm=T), length=100)
	# STILL WITH scale()
Turbid <- c(0,1)
# GS <- scale(c(2,6,10))
Scaled_group_sizes <- sort(unique(scale(Fdata$Group_size)))
GS <- Scaled_group_sizes

Plot_data <- expand.grid(Behaviour,Turbid,GS)
colnames(Plot_data) <- c("Behaviour","Turbid","GS")

Predicted <- fixef(m21)[9]*Plot_data$Behaviour*Plot_data$GS +
	fixef(m21)[8]*Plot_data$Behaviour*Plot_data$Turbid +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*Plot_data$GS +
	fixef(m21)[3]*Plot_data$Turbid +
	fixef(m21)[2]*Plot_data$Behaviour +
	fixef(m21)[1]	
range(predict(m21))
range(Predicted)
	# CHECK: CLOSE MATCH
range(log(Fdata$Latency_of_attack))

Plot_data$Behaviour_not_scaled <- rep( seq(from = min(Fdata$Med_group_polarisation,na.rm=T), to = max(Fdata$Med_group_polarisation,na.rm=T), length=100) ,6)
	# FOR PLOTTING ONLY
# plot(Plot_data$Behaviour_not_scaled~Plot_data$Behaviour)
	# CHECK: IS A LINEAR SCALING

par(mar=c(3, 1, 1.3, 0.4))
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]~Fdata$Med_group_polarisation[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=c(0,1)
	,xlab="",ylab="",main="2 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==2]~Fdata$Med_group_polarisation[Fdata$Treatment=="Clear" & Fdata$Group_size==2]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[1]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[1]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[1]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[1]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,0.5,1),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("C",side=3,line=0.07,cex=1,adj=0,font=2)

par(mar=c(3, 1, 1.3, 0.4))		
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]~Fdata$Med_group_polarisation[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=c(0,1)
	,xlab="",ylab="",main="6 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==6]~Fdata$Med_group_polarisation[Fdata$Treatment=="Clear" & Fdata$Group_size==6]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[2]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[2]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[2]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[2]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,0.5,1),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
par("cex")
mtext("Median group polarisation",side=1,line=1.7,adj=0.5,outer=F,cex=0.8)

par(mar=c(3, 1, 1.3, 0.4))
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==10]~Fdata$Med_group_polarisation[Fdata$Treatment=="Turbid" & Fdata$Group_size==10]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=c(0,1)
	,xlab="",ylab="",main="10 fish",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==10]~Fdata$Med_group_polarisation[Fdata$Treatment=="Clear" & Fdata$Group_size==10]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[3]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==0 & Plot_data$GS==Scaled_group_sizes[3]],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[3]])~Plot_data$Behaviour_not_scaled[Plot_data$Turbid==1 & Plot_data$GS==Scaled_group_sizes[3]],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,0.5,1),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)

dev.off()


# Med_speed
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Med_speed)*Treatment+scale(Med_speed)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Med_speed)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Med_speed)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Med_speed)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Med_speed)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3 WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY IMPORTANT INTERACTION
    dAICc df
m23  0.0  11
m21  0.1  12
m25  2.4  14
m24  3.7  10
m22  3.9  11
m2  53.2  9 
length(residuals(m2))


# Med_bearing_of_NN
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Med_bearing_of_NN)*Treatment+scale(Med_bearing_of_NN)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Med_bearing_of_NN)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Med_bearing_of_NN)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Med_bearing_of_NN)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Med_bearing_of_NN)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3 WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# WITHOUT Med_bearing_of_NN IS BEST MODEL
	# NO INTERACTIONS IMPROVE THE MODEL FIT
    dAICc df
m22  0.0  11
m2   0.9  9 
m21  1.5  12
m24  1.8  10
m23  3.4  11
m25  4.0  14
length(residuals(m2))


# Min_dist_stimulus
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)*Treatment+scale(Min_dist_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Min_dist_stimulus)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Min_dist_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Min_dist_stimulus)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Min_dist_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3 WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
# plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# MODEL WITH BEHAVIOURAL VARIABLE x TURBIDITY IMPROVES FIT
	# INCLUDING BEHAVIOURAL VARIABLE x GROUP SIZE WORSENS MODEL FIT
    dAICc df
m22  0.0  11
m21  1.0  12
m25  5.0  14
m24  9.8  10
m23 10.7  11
m2  23.4  9 
length(residuals(m2))

summary(m22)
confint(m22)

fixef(m22)
Behaviour <- rep( seq(from = min(scale(Fdata$Min_dist_stimulus_cm),na.rm=T), to = max(scale(Fdata$Min_dist_stimulus_cm),na.rm=T), length=100) ,2)
	# STILL WITH scale()
Turbid <- c(rep(0,100),rep(1,100))
Predicted <- fixef(m22)[8]*Behaviour*Turbid +
	fixef(m22)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m22)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m22)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m22)[4]*mean(scale(Fdata$Group_size)) +
	fixef(m22)[3]*Turbid +
	fixef(m22)[2]*Behaviour +
	fixef(m22)[1]
range(predict(m22))
range(Predicted)

Behaviour_not_scaled <- rep( seq(from = min(Fdata$Min_dist_stimulus_cm,na.rm=T), to = max(Fdata$Min_dist_stimulus_cm,na.rm=T), length=100) ,2)
	# FOR PLOTTING ONLY
# plot(Behaviour_not_scaled~Behaviour)
	# CHECK: IS A LINEAR SCALING

# FIG. 4

tiff("Fig 4.tiff", width=8.7,height=11, units = "cm", type="windows", compression="lzw", res=600)
par(oma=c(0, 3, 0, 0))
Panels <- layout( matrix(c(1,2,3,3,4,4),nrow=3,byrow=T),heights = c(1/3,1/3,1/3) )
	# THIRD ROW ISN'T USED, BUT LATER CROPPED IN ANOTHER SOFTWARE. ENSURES THAT PLOT HEIGHTS ARE CONSISTENT BETWEEN FIGURES.
# layout.show( Panels )
Panels

par(mar=c(1, 1, 1.3, 0.4))
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]~Fdata$Min_dist_stimulus_cm[Fdata$Treatment=="Turbid" & Fdata$Group_size==2]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T),xlim=c(0,max(Fdata$Min_dist_stimulus_cm,na.rm=T))
	,xlab="",ylab="",main="",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==2]~Fdata$Min_dist_stimulus_cm[Fdata$Treatment=="Clear" & Fdata$Group_size==2]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Turbid==0])~Behaviour_not_scaled[Turbid==0],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Turbid==1])~Behaviour_not_scaled[Turbid==1],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,40,80,120),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("A",side=3,line=0.07,adj=0,font=2)
mtext("Minimum distance\nto stimulus (cm)",side=1,line=2.7,adj=0.5,outer=F,cex=0.8)
mtext("Latency to attack (s)",side=2,line=2.3,adj=0.5,outer=F,cex=0.8)

	
# Min_bearing_to_stimulus
m2<-lmer(log(Latency_of_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# AS m2 ABOVE
m21<-lmer(log(Latency_of_attack)~scale(Min_bearing_to_stimulus)*Treatment+scale(Min_bearing_to_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BOTH INTERACTIONS
m22<-lmer(log(Latency_of_attack)~scale(Min_bearing_to_stimulus)*Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x TURBIDITY ONLY
m23<-lmer(log(Latency_of_attack)~Treatment+scale(Min_bearing_to_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# BEHAVIOURAL VARIABLE x GROUP SIZE ONLY
m24<-lmer(log(Latency_of_attack)~Treatment+scale(Min_bearing_to_stimulus)+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# MAIN EFFECTS ONLY
m25<-lmer(log(Latency_of_attack)~Treatment*scale(Min_bearing_to_stimulus)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
	# 3 WAY INTERACTION
simulationOutput <- simulateResiduals(m25,n=1000)
# plot(simulationOutput)
ICtab(m2,m21,m22,m23,m24,m25,type=c("AICc"))
	# MODEL WITH BEHAVIOURAL VARIABLE x TURBIDITY IMPROVES FIT; BOTH TWO-WAY INTERACTIONS MODEL IS NOT BETTER ENOUGH (ONLY 0.1 UNITS)
    dAICc df
m22  0.0  11
m21  1.5  12
m24  4.1  10
m25  4.7  14
m23  5.8  11
m2  47.2  9 
length(residuals(m2))

summary(m22)
confint(m22)

fixef(m22)
Behaviour <- rep( seq(from = min(scale(Fdata$Min_bearing_to_stimulus),na.rm=T), to = max(scale(Fdata$Min_bearing_to_stimulus),na.rm=T), length=100) ,2)
	# STILL WITH scale()
Turbid <- c(rep(0,100),rep(1,100))
Predicted <- fixef(m22)[8]*Behaviour*Turbid +
	fixef(m22)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m22)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m22)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m22)[4]*mean(scale(Fdata$Group_size)) +
	fixef(m22)[3]*Turbid +
	fixef(m22)[2]*Behaviour +
	fixef(m22)[1]
range(predict(m22))
range(Predicted)

Behaviour_not_scaled <- rep( seq(from = min(Fdata$Min_bearing_to_stimulus,na.rm=T), to = max(Fdata$Min_bearing_to_stimulus,na.rm=T), length=100) ,2)
	# FOR PLOTTING ONLY
# plot(Behaviour_not_scaled~Behaviour)
	# CHECK: IS A LINEAR SCALING

par(mar=c(1, 1, 1.3, 0.4))		
plot(Fdata$Latency_of_attack[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]~Fdata$Min_bearing_to_stimulus[Fdata$Treatment=="Turbid" & Fdata$Group_size==6]
	,col=rgb(128/255,128/255,128/255,alpha=0.4),pch=16,cex=0.8
	,log="y",xaxt="n",yaxt="n"
	,ylim=range(Fdata$Latency_of_attack,na.rm=T)
	,xlab="",ylab="",main="",cex.main=1,font.main=1,las=1)
points(Fdata$Latency_of_attack[Fdata$Treatment=="Clear" & Fdata$Group_size==6]~Fdata$Min_bearing_to_stimulus[Fdata$Treatment=="Clear" & Fdata$Group_size==6]
	,col=rgb(0/255,204/255,204/255,alpha=1),cex=0.3)
lines(exp(Predicted[Turbid==0])~Behaviour_not_scaled[Turbid==0],col=rgb(0/255,204/255,204/255),lwd=1)
lines(exp(Predicted[Turbid==1])~Behaviour_not_scaled[Turbid==1],col=rgb(128/255,128/255,128/255),lwd=1)
axis(side=1,at=c(0,40,80,120),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("B",side=3,line=0.07,adj=0,font=2)
mtext("Minimum bearing to\nstimulus (degrees)",side=1,line=2.7,adj=0.5,outer=F,cex=0.8)

par(mar=c(0.3, 1, 4, 0.4))
boxplot(Fdata$Tortuosity_of_first_to_attack~Fdata$Treatment*Fdata$Group_size,log="y",las=1,xlab="",ylab="",xaxt="n",yaxt="n"
	,col=rep(c(rgb(0/255,204/255,204/255,alpha=0),rgb(128/255,128/255,128/255,alpha=0.4)),3)
	,border=rep(c(rgb(0/255,204/255,204/255),"black"),3))
axis(side=1,at=1:6,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=1,at=c(1.5,3.5,5.5),labels=F,mgp=c(1.5,0.6,0),las=1,tcl=0,line=1,lty="blank")
axis(side=2,at=c(2,5,20,50,200),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("C",side=3,line=0.07,adj=0,font=2)
mtext("Tortuosity of\nfirst to attack",side=2,line=1.8,adj=0.5,outer=F,cex=0.8)

par(mar=c(2.6, 1, 1.7, 0.4))
boxplot(Fdata$Med_speed_of_first_to_reach_cm~Fdata$Treatment*Fdata$Group_size,las=1,xlab="",ylab="",xaxt="n",yaxt="n"
	,col=rep(c(rgb(0/255,204/255,204/255,alpha=0),rgb(128/255,128/255,128/255,alpha=0.4)),3)
	,border=rep(c(rgb(0/255,204/255,204/255),"black"),3))		
axis(side=1,at=1:6,labels=rep(c("Clear","Turbid"),3),mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=1,at=c(1.5,3.5,5.5),labels=c("2 fish","6 fish","10 fish"),mgp=c(1.5,0.6,0),las=1,tcl=0,line=1,lty="blank")
axis(side=2,at=c(0,10,20,30),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("D",side=3,line=0.07,adj=0,font=2)
mtext("Median speed of\nfirst to attack (cm/s)",side=2,line=1.8,adj=0.5,outer=F,cex=0.8)

dev.off()


## DOES TURBIDITY AND GROUP SIZE EFFECT BEHAVIOUR IMMEDIATELY BEFORE THE PRESENTATION APPEARS?

# Med_distance_group_centroid
hist(Fdata$Med_distance_group_centroid_cm)

m1<-glmer.nb(Med_distance_group_centroid_cm~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
simulationOutput <- simulateResiduals(m1,n=1000)
testDispersion(simulationOutput)
	# DISPERSION OK WITH NEG BIN
m2<-glmer.nb(Med_distance_group_centroid_cm~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
m3<-glmer.nb(Med_distance_group_centroid_cm~Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
m4<-glmer.nb(Med_distance_group_centroid_cm~          scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
m5<-glmer.nb(Med_distance_group_centroid_cm~          			   scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
ICtab(m1,m2,m3,m4,m5,type=c("AICc"))
	# MAIN EFFECTS ONLY IS THE BEST SUPPORTED MODEL
   dAICc df
m1  0.0  10
m2  1.8  9 
m4 26.8  8 
m3 52.3  8 
m5 69.0  7
length(residuals(m2))

summary(m2)
confint.merMod(m2,method = "Wald")


# Med_group_polarisation
hist(Fdata$Med_group_polarisation)

Fdata$Polarisation_trans <- 1-Fdata$Med_group_polarisation
hist(Fdata$Polarisation_trans)
m1<-lmer(Med_group_polarisation~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
simulationOutput <- simulateResiduals(m1,n=1000)
plot(simulationOutput)
	# ASSUMPTIONS NOT MET
m1<-glmer.nb(Polarisation_trans~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
simulationOutput <- simulateResiduals(m1,n=1000)
testDispersion(simulationOutput)
	# NEG BIN ASSUMPTION NOT MET
m1<-glmer(Med_group_polarisation~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata, family="poisson")
simulationOutput <- simulateResiduals(m1,n=1000)
testDispersion(simulationOutput)
	# POISSON ASSUMPTION NOT MET

range(Fdata$Med_group_polarisation)
	# DOESN'T NEED TRANSFORMING FOR BETA REGRESSION
m1<-glmmTMB(Med_group_polarisation~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
simulationOutput <- simulateResiduals(m1,n=1000)
plot(simulationOutput)
	# ASSUPTIONS GOOD WITH BETA
m2<-glmmTMB(Med_group_polarisation~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
m3<-glmmTMB(Med_group_polarisation~Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
m4<-glmmTMB(Med_group_polarisation~          scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
m5<-glmmTMB(Med_group_polarisation~                            scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
ICtab(m1,m2,m3,m4,m5,type=c("AICc"))
	# MAIN EFFECTS ONLY MODEL IS THE MOST LIKELY MODEL
   dAICc df
m2  0.0  9 
m1  1.5  10
m4 11.1  8 
m3 67.3  8 
m5 72.4  7 
length(residuals(m2))
summary(m2)
	# NEGATIVE EFFECT OF GROUP SIZE; LESS TORTUOUS IN LARGER GROUPS
confint(m2)
	# 95% CONFIDENCE INTERVALS FOR COEFFICIENT ESTIMATES
	
# FIG. 3

tiff("Fig 3.tiff", width=8.7,height=11, units = "cm", type="windows", compression="lzw", res=600)
par(oma=c(0, 3, 0, 0))
par(mar=c(2, 1, 2.3, 0.4))
layout( matrix(c(1,2,3),nrow=3,byrow=T) )
	# THIRD ROW ISN'T USED, BUT LATER CROPPED IN ANOTHER SOFTWARE. ENSURES THAT PLOT HEIGHTS ARE CONSISTENT BETWEEN FIGURES.

boxplot(Fdata$Med_distance_group_centroid_cm~Fdata$Treatment*Fdata$Group_size,las=1,xlab="",ylab="",xaxt="n",yaxt="n"
	,col=rep(c(rgb(0/255,204/255,204/255,alpha=0),rgb(128/255,128/255,128/255,alpha=0.4)),3)
	,border=rep(c(rgb(0/255,204/255,204/255),"black"),3))
axis(side=1,at=1:6,labels=F,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=2,at=c(0,20,40,60),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("A",side=3,line=0.07,cex=1,adj=0,font=2)
mtext("Median distance to\ngroup centroid (cm)",side=2,line=1.8,adj=0.5,outer=F,cex=0.8)

par(mar=c(4, 1, 0.3, 0.4))
boxplot(Fdata$Med_group_polarisation~Fdata$Treatment*Fdata$Group_size,las=1,xlab="",ylab="",xaxt="n",yaxt="n",ylim=c(0,1)
	,col=rep(c(rgb(0/255,204/255,204/255,alpha=0),rgb(128/255,128/255,128/255,alpha=0.4)),3)
	,border=rep(c(rgb(0/255,204/255,204/255),"black"),3))
axis(side=1,at=1:6,labels=rep(c("Clear","Turbid"),3),mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
axis(side=1,at=c(1.5,3.5,5.5),labels=c("2 fish","6 fish","10 fish"),mgp=c(1.5,0.6,0),las=1,tcl=0,line=1,lty="blank")
axis(side=2,at=c(0,0.5,1),labels=T,mgp=c(1.5,0.6,0),las=1,tcl=-0.4)
mtext("B",side=3,line=0.07,cex=1,adj=0,font=2)
mtext("Median group\npolarisation",side=2,line=1.8,adj=0.5,outer=F,cex=0.8)

dev.off()

plot(Fdata$Med_group_polarisation ~ Fdata$Med_distance_group_centroid_cm)


## WHAT IS THE PREDICTED IMPACT OF THESE CHANGES OF COLLECTIVE BEHAVIOUR ON THE LATENCY TO ATTACK?

hist(Fdata$Med_distance_group_centroid_cm)
m2<-glmer.nb(Med_distance_group_centroid_cm~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
fixef(m2)
	# ALL FIXED EFFECT COEFFICIENTS
range(predict(m2,type="response"))
predicted_distance_centroid_clear <- exp( fixef(m2)[1] +
	fixef(m2)[2]*0 +
	fixef(m2)[3]*sort(unique(scale(Fdata$Group_size))) +
		# 1ST IS GROUP SIZE 2, 2ND SIZE 6, 3RD SIZE 10
	fixef(m2)[4]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m2)[5]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m2)[6]*mean(scale(Fdata$Day_to_save)) )
predicted_distance_centroid_turbid <- exp( fixef(m2)[1] +
	fixef(m2)[2]*1 +
	fixef(m2)[3]*sort(unique(scale(Fdata$Group_size))) +
		# 1ST IS GROUP SIZE 2, 2ND SIZE 6, 3RD SIZE 10
	fixef(m2)[4]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m2)[5]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m2)[6]*mean(scale(Fdata$Day_to_save)) )
		# 0 OR 1 FOR COEFFICIENTS RELATED TO TURBIDITY, AS CLEAR = 0 AND TURBID = 1

plot(scale(Fdata$Med_distance_group_centroid_cm)~Fdata$Med_distance_group_centroid_cm)
predicted_distance_centroid_clear_scaled  <- (predicted_distance_centroid_clear  - mean(Fdata$Med_distance_group_centroid_cm,na.rm=T)) / sd(Fdata$Med_distance_group_centroid_cm,na.rm=T)
predicted_distance_centroid_turbid_scaled <- (predicted_distance_centroid_turbid - mean(Fdata$Med_distance_group_centroid_cm,na.rm=T)) / sd(Fdata$Med_distance_group_centroid_cm,na.rm=T)
	# CONVERTING DISTANCE TO GROUP CENTROID ON ORIGINAL SCALE TO THE SCALE USED IN THE LINEAR MODEL
points(predicted_distance_centroid_clear_scaled~predicted_distance_centroid_clear,col="red",pch=16)
points(predicted_distance_centroid_turbid_scaled~predicted_distance_centroid_turbid,col="green",pch=16)
	# CHECK: SHOULD FALL ON LINE MADE BY POINTS

m21<-lmer(log(Latency_of_attack)~scale(Med_distance_group_centroid_cm)*Treatment+scale(Med_distance_group_centroid_cm)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
predicted_latency_clear <- exp( fixef(m21)[9]*predicted_distance_centroid_clear_scaled*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[8]*predicted_distance_centroid_clear_scaled*1 +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[3]*1 +
	fixef(m21)[2]*predicted_distance_centroid_clear_scaled +
	fixef(m21)[1] )
		# THE LATENCY IN TURBID WATER IF THE FISH HAD THE COHESION THEY WOULD HAVE IN CLEAR WATER (AT EACH GROUP SIZE) 
predicted_latency_turbid <- exp( fixef(m21)[9]*predicted_distance_centroid_turbid_scaled*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[8]*predicted_distance_centroid_turbid_scaled*1 +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[3]*1 +
	fixef(m21)[2]*predicted_distance_centroid_turbid_scaled +
	fixef(m21)[1] )
data.frame(predicted_distance_centroid_clear,predicted_distance_centroid_turbid,predicted_distance_centroid_clear_scaled,predicted_distance_centroid_turbid_scaled,predicted_latency_clear,predicted_latency_turbid)
(predicted_latency_clear-predicted_latency_turbid)/predicted_latency_clear

Latency_improvement <- predicted_latency_clear-predicted_latency_turbid
data.frame(predicted_distance_centroid_clear,predicted_distance_centroid_turbid
	,predicted_latency_clear,predicted_latency_turbid
	,Latency_improvement
	,Latency_improvement*100/predicted_latency_clear)


m2<-glmmTMB(Med_group_polarisation~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=beta_family(link="logit"))
fixef(m2)$cond
	# ALL FIXED EFFECT COEFFICIENTS
range(predict(m2,type="response"))
predicted_polarisation_clear <- inv.logit( fixef(m2)$cond[1] +
	fixef(m2)$cond[2]*0 +
	fixef(m2)$cond[3]*sort(unique(scale(Fdata$Group_size))) +
		# 1ST IS GROUP SIZE 2, 2ND SIZE 6, 3RD SIZE 10
	fixef(m2)$cond[4]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m2)$cond[5]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m2)$cond[6]*mean(scale(Fdata$Day_to_save)) )
predicted_polarisation_turbid <- inv.logit( fixef(m2)$cond[1] +
	fixef(m2)$cond[2]*1 +
	fixef(m2)$cond[3]*sort(unique(scale(Fdata$Group_size))) +
		# 1ST IS GROUP SIZE 2, 2ND SIZE 6, 3RD SIZE 10
	fixef(m2)$cond[4]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m2)$cond[5]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m2)$cond[6]*mean(scale(Fdata$Day_to_save)) )
		# 0 OR 1 FOR COEFFICIENTS RELATED TO TURBIDITY, AS CLEAR = 0 AND TURBID = 1

plot(scale(Fdata$Med_group_polarisation)~Fdata$Med_group_polarisation)
predicted_polarisation_clear_scaled  <- (predicted_polarisation_clear  - mean(Fdata$Med_group_polarisation,na.rm=T)) / sd(Fdata$Med_group_polarisation,na.rm=T)
predicted_polarisation_turbid_scaled <- (predicted_polarisation_turbid - mean(Fdata$Med_group_polarisation,na.rm=T)) / sd(Fdata$Med_group_polarisation,na.rm=T)
	# CONVERTING POLARISATION ON ORIGINAL SCALE TO THE SCALE USED IN THE LINEAR MODEL
points(predicted_polarisation_clear_scaled~predicted_polarisation_clear,col="red",pch=16)
points(predicted_polarisation_turbid_scaled~predicted_polarisation_turbid,col="green",pch=16)
	# CHECK - SHOULD FALL ON LINE MADE BY POINTS

m21<-lmer(log(Latency_of_attack)~scale(Med_group_polarisation)*Treatment+scale(Med_group_polarisation)*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata,REML=F)
predicted_latency_clear <- exp( fixef(m21)[9]*predicted_polarisation_clear_scaled*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[8]*predicted_polarisation_clear_scaled*1 +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[3]*1 +
	fixef(m21)[2]*predicted_polarisation_clear_scaled +
	fixef(m21)[1] )
		# THE LATENCY IN TURBID WATER IF THE FISH HAD THE POLARISATION THEY WOULD HAVE IN CLEAR WATER (AT EACH GROUP SIZE) 
predicted_latency_turbid <- exp( fixef(m21)[9]*predicted_polarisation_turbid_scaled*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[8]*predicted_polarisation_turbid_scaled*1 +
	fixef(m21)[7]*mean(scale(Fdata$Day_to_save)) +
	fixef(m21)[6]*mean(scale(Fdata$Trial_in_day)) +
	fixef(m21)[5]*mean(scale(Fdata$Presentation_in_trial)) +
	fixef(m21)[4]*sort(unique(scale(Fdata$Group_size))) +
	fixef(m21)[3]*1 +
	fixef(m21)[2]*predicted_polarisation_turbid_scaled +
	fixef(m21)[1] )

Latency_improvement <- predicted_latency_clear-predicted_latency_turbid
data.frame(predicted_polarisation_clear,predicted_polarisation_turbid
	# ,predicted_polarisation_clear_scaled,predicted_polarisation_turbid_scaled
	,predicted_latency_clear,predicted_latency_turbid
	,Latency_improvement
	,Latency_improvement*100/predicted_latency_clear)


## DOES TURBIDITY AND GROUP SIZE EFFECT TORTUOSITY AFTER THE PRESENTATION AND BEFORE THE STIMULUS IS REACHED?

str(Fdata)

# Tortuosity_of_first_to_attack
hist(Fdata$Tortuosity_of_first_to_attack)
hist(log(Fdata$Tortuosity_of_first_to_attack))
range(Fdata$Tortuosity_of_first_to_attack)

m1<-glmer.nb(Tortuosity_of_first_to_attack~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number), data=Fdata)
simulationOutput <- simulateResiduals(m1,n=1000)
testDispersion(simulationOutput)
	# MASSIVELY OVER DISPERSED

m1<-glmmTMB(log(Tortuosity_of_first_to_attack)~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=Gamma(link="log"))
simulationOutput <- simulateResiduals(m1,n=1000)
plot(simulationOutput)
	# WHEN LOGGED, DIAGNOSTICS LOOK GOOD WITH GAMMA
m2<-glmmTMB(log(Tortuosity_of_first_to_attack)~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=Gamma(link="log"))
m3<-glmmTMB(log(Tortuosity_of_first_to_attack)~Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=Gamma(link="log"))
m4<-glmmTMB(log(Tortuosity_of_first_to_attack)~          scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=Gamma(link="log"))
m5<-glmmTMB(log(Tortuosity_of_first_to_attack)~                            scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,family=Gamma(link="log"))
ICtab(m1,m2,m3,m4,m5,type=c("AICc"))
	# MAIN EFFECTS ONLY IS THE BEST SUPPORTED MODEL
   dAICc df
m2  0.0  9 
m1  0.1  10
m4 29.1  8 
m3 52.4  8 
m5 71.7  7 
length(residuals(m2))

summary(m2)
confint(m2)
	

# Med_speed_of_first_to_reach
hist(Fdata$Med_speed_of_first_to_reach)
range(Fdata$Med_speed_of_first_to_reach)

m1<-lmer(Med_speed_of_first_to_reach~Treatment*scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,REML=F)
simulationOutput <- simulateResiduals(m1,n=1000)
plot(simulationOutput)
m2<-lmer(Med_speed_of_first_to_reach~Treatment+scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,REML=F)
m3<-lmer(Med_speed_of_first_to_reach~Treatment                  +scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,REML=F)
m4<-lmer(Med_speed_of_first_to_reach~          scale(Group_size)+scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,REML=F)
m5<-lmer(Med_speed_of_first_to_reach~                            scale(Presentation_in_trial)+scale(Trial_in_day)+scale(Day_to_save)+(1|Tank/Trial_number),data=Fdata,REML=F)
ICtab(m1,m2,m3,m4,m5,type=c("AICc"))
   dAICc df
m1  0.0  10
m2  6.3  9 
m4 16.7  8 
m3 47.2  8 
m5 53.9  7 
length(residuals(m2))

summary(m1)
confint(m1)

plot(log(Fdata$Tortuosity_of_first_to_attack)~Fdata$Med_speed_of_first_to_reach)