---
title: Kinematic performance declines as group size increases during escape responses
  in a schooling coral reef fish
author: "Ms Monica D Bacchus & Dr Lauren E Nadler"
date: "2023-09-12"
output: html_document
---
Escaping predation is essential for species survival, but prey must effectively match their response to the perceived threat imposed by a predator. For social animals, one mechanism to reduce risk of predation is living in larger group sizes, which dilutes each individual’s risk of capture. When a predator attacks, individuals from a range of taxa (e.g., fishes, sharks, amphibians) use the rapid, anaerobically-fueled burst swimming behavior, known as the fast-start response, to evade the attack. Here, using the schooling coral reef damselfish Chromis viridis, we assess if there is an optimal group size that maximizes both individual escape response as well as group cohesion and coordination following a simulated predator attack, comparing schools composed of four, eight, and sixteen fish. We found that fish in various group sizes exhibited no difference in their reaction timing to a simulated predator attack (i.e., latency), but larger groups exhibited a slower kinematic response (i.e., lower average turning rate and shorter distance covered during the escape response), potentially because larger groups perceived the predator attack as less risky due to safety in numbers. Both school cohesion and coordination (as measured through alignment and nearest neighbor distance, respectively) declined in the 100ms after the predator’s attack. While there was no impact of group size on alignment, larger group sizes exhibited closer nearest neighbor distance at all time stamps. This study highlights that larger group sizes may allow individuals to save energy on costly behavioral responses to avoid predators, potentially through a greater threshold of threat necessary to trigger a rapid escape response. 

```{r, message = F}
#LIBRARIES
library(lme4)
library(MuMIn)
library(emmeans)
library(car)
library(ggplot2)
```

# LATENCY

Latency, or reaction timing, was measured as the time interval (in milliseconds) between the aerial mechanical stimulus first breaking the water's surface and the fish's initial head movement. Non-responders (n=3, 0.8% of all fish included in this study; those fish that did not respond within two seconds of the stimulus) were assigned the maximum measured latency in this study (1003.8 ms). 

```{r, message = F, warning=FALSE}
fs_all <- read.csv("/Users/laurennadler/Desktop/MyRData/Monica_groupsize_faststart.csv")
fs_all$Group.Size = as.factor(fs_all$Group.Size)
lat.hist <- ggplot(fs_all, aes(x=Latency, color = Group.Size, fill = Group.Size)) + 
  geom_histogram(position = "dodge", alpha = 0.5, bins = 6) +
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.title=element_blank()) +
  theme(legend.position="top") 
lat.hist

#latency model
allLat.lmer = lmer((Latency + 0.001)~ Group.Size*Stimulus.distance + (1|Video), data=fs_all, na.action = "na.omit") 
plot(allLat.lmer) #highly heteroscedastic
qqnorm(resid(allLat.lmer)) 
qqline(resid(allLat.lmer))
shapiro.test(resid(allLat.lmer)) #p=9.752e-09 so non normal
bartlett.test(resid(allLat.lmer), fs_all$Group.Size) #p=0.03 so heteroscedastic


allLat.lmer.bc <- powerTransform(allLat.lmer, family="bcPower") #used boxcox transformation due to violation of assumptions of normality and homoscedasticity
summary(allLat.lmer.bc) #both significant
allLat.lmer.bc$roundlam


allLat.bc.lmer = lmer(bcPower((Latency+0.001), allLat.lmer.bc$roundlam) ~  Group.Size*Stimulus.distance + (1|Video), data=fs_all, na.action = "na.omit")
plot(allLat.bc.lmer)
qqnorm(resid(allLat.bc.lmer)) 
qqline(resid(allLat.bc.lmer))
shapiro.test(resid(allLat.bc.lmer)) #still p < 0.05 but visual inspection of residuals plot dramatically improved by boxcox transformation so continuing with lmer with boxcox transformation
bartlett.test(resid(allLat.bc.lmer), fs_all$Group.Size) 

allLat.bc.lmer2 = lmer(bcPower((Latency+0.001), allLat.lmer.bc$roundlam) ~  Group.Size+Stimulus.distance + (1|Video), data=fs_all, na.action = "na.omit")
anova(allLat.bc.lmer,allLat.bc.lmer2) #more complex model has a lower AIC so maintaining allLat.bc.lmer

summary(allLat.bc.lmer)
Anova(allLat.bc.lmer, test="F")

r.squaredGLMM(allLat.bc.lmer)
latency.plot <- ggplot(fs_all, aes(x=Stimulus.distance, y=Latency, colour=Group.Size)) + 
  geom_point(aes(color = Group.Size), alpha = 0.7, size=1) + 
  geom_smooth(method="glm", fullrange=TRUE, aes(fill = Group.Size)) + 
  xlab(bquote('Stimulus distance (mm)')) + ylab(bquote('Latency (ms)')) + 
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.title=element_blank()) +
  scale_y_log10() +
  theme(legend.position = "top")
latency.plot
```

# AVERAGE TURNING RATE & DISTANCE COVERED 

We quantified fast-start characteristics associated with the initial unilateral body bend following stimulation (i.e., stage 1 of the fast-start) and the subsequent contralateral body bend (i.e., stage 2 of the fast-start). Individual fast-start kinematic performance was characterized through average turning rate (AVT, the maximum turning angle achieved by the fish during stage 1 divided by the time it took to achieve that angle, which serves as a proxy for the response’s agility through speed of muscle contraction), and distance covered (DC; distance moved in the first 42 ms of the reaction, which is the average time for this species to achieve stages 1 and 2; used as a proxy for swimming speed and acceleration). Since these traits are influenced by the stimulus distance (the distance between the fish’s center of mass and the stimulus), this trait was also measured and included as a covariate in all analyses (including latency above). 

```{r, message = F, warning=FALSE}
#avt - NAs removed (ie - fish too close to wall and non-responders) 
fs_noNA <- read.csv("/Users/laurennadler/Desktop/MyRData/Monica_groupsize_faststart no nas.csv")
fs_noNA$Group.Size = as.factor(fs_noNA$Group.Size)
avt.lmer = lmer(Average.Turning.Rate.per.s ~ Group.Size*Stimulus.distance + (1|Video), data=fs_noNA, na.action = "na.omit")
plot(avt.lmer) 
qqnorm(resid(avt.lmer))
qqline(resid(avt.lmer))
shapiro.test(resid(avt.lmer)) #p = 0.003 suggesting non normal but disregarding because the qqplot looks fine
bartlett.test(resid(avt.lmer), fs_noNA$Group.Size)
avt.lmer2 = lmer(Average.Turning.Rate.per.s ~ Group.Size+Stimulus.distance + (1|Video), data=fs_noNA, na.action = "na.omit")
anova(avt.lmer,avt.lmer2) #suggests that the more complex lmer better
summary(avt.lmer)
Anova(avt.lmer, test="F") 
r.squaredGLMM(avt.lmer)
emmeans(avt.lmer, pairwise~Group.Size, adjust = c("Tukey")) 
avt.plot <- ggplot(fs_noNA, aes(x=Stimulus.distance, y=Average.Turning.Rate.per.s, colour=Group.Size)) + 
  geom_point(aes(color = Group.Size), alpha = 0.7, size=1) + 
  geom_smooth(method="glm", fullrange=TRUE, aes(fill = Group.Size)) + 
  xlab(bquote('Stimulus distance (cm)')) + ylab(bquote('ATR (o/s)')) + 
  theme_classic(base_family='Arial', base_size = 20) +
  theme(legend.title=element_blank()) +
  theme(legend.position = "top") 
avt.plot
avt.emm <-emmeans(avt.lmer, ~Group.Size) 
summary(avt.emm)
avt.emm<-as.data.frame(avt.emm)
avt.emplot = ggplot(avt.emm, aes(x=Group.Size, y=emmean, color=Group.Size, fill=Group.Size)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL, color = Group.Size), na.rm = TRUE, shape=21, size=.7, data = avt.emm) + 
  geom_jitter(aes(x=Group.Size, y = Average.Turning.Rate.per.s, group=Group.Size), shape = 21, size=1.2, width=.25, alpha = 0.35, data = fs_noNA) + 
  xlab(bquote('Group Size')) + ylab(bquote('ATR (o/s)')) + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +  
  theme(axis.text = element_text(size = 12)) + theme(axis.title = element_text(size = 12)) +
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.position = "none") 

avt.emplot

#DC 
dc.lmer = lmer(Distance.covered ~ Group.Size*Stimulus.distance + (1|Video), data=fs_noNA, na.action = "na.omit")
plot(dc.lmer)
qqnorm(resid(dc.lmer)) #needs a boxcox transformation
qqline(resid(dc.lmer))
shapiro.test(resid(dc.lmer)) #p-value = 1.902e-08 so non normal
bartlett.test(resid(dc.lmer), fs_noNA$Group.Size) # p-value = 3.129e-05 so heteroscedastic
dc.lmer.bc <- powerTransform(dc.lmer, family="bcPower") 
summary(dc.lmer.bc) #both significant

dc.lmer.bc$roundlam

dc.bc.lmer = lmer(bcPower(Distance.covered, dc.lmer.bc$roundlam) ~  Group.Size*Stimulus.distance + (1|Video), data=fs_noNA, na.action = "na.omit")
plot(dc.bc.lmer)
qqnorm(resid(dc.bc.lmer)) 
qqline(resid(dc.bc.lmer))
shapiro.test(resid(dc.bc.lmer))  #p=0.0001244 - glmer robust to minor deviations in assumptions so continuing with this transformed/improved model
bartlett.test(resid(dc.bc.lmer), fs_noNA$Group.Size) #p= 0.002 but visual inspection of the plot looks fine so disregarding

dc.bc.lmer2 = lmer(bcPower(Distance.covered, dc.lmer.bc$roundlam) ~  Group.Size+Stimulus.distance + (1|Video), data=fs_noNA, na.action = "na.omit")
anova(dc.bc.lmer, dc.bc.lmer2) #simpler lmer2 only has an AIC that is lower by <1 so sticking with most complex model
summary(dc.bc.lmer)
Anova(dc.bc.lmer, test="F")
r.squaredGLMM(dc.bc.lmer)
backtrans.dc <- update(ref_grid(dc.bc.lmer), tran = make.tran("boxcox", 0.5)) #back transforms from boxcox so emmeans produces sensible estimates of means and SE where 0.5 is the lambda estimated by powertransform function above
emmeans(backtrans.dc, pairwise~Group.Size, adjust = c("Tukey"), type = "response") 
dc.emm<-emmeans(backtrans.dc, ~Group.Size, type = "response") 
summary(dc.emm)
dc.emm<-as.data.frame(dc.emm)
dc.emplot = ggplot(dc.emm, aes(x=Group.Size, y=response, color=Group.Size, fill=Group.Size)) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL, color = Group.Size), na.rm = TRUE, shape=21, size=.7, data = dc.emm) + 
  geom_jitter(aes(x=Group.Size, y = Distance.covered, group=Group.Size), shape = 21, size=1.2, width=.25, alpha = 0.35, data = fs_noNA) + 
  xlab(bquote('Group Size')) + ylab(bquote('Distance covered (mm)')) + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +  
  theme(axis.text = element_text(size = 12)) + theme(axis.title = element_text(size = 12)) +
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.position = "none") 

dc.emplot

dc.plot <- ggplot(fs_noNA, aes(x=Stimulus.distance, y=Distance.covered, colour=Group.Size)) + 
  geom_point(aes(color = Group.Size), alpha = 0.7, shape=21, size=.7) + 
  geom_smooth(method="glm", fullrange=TRUE, aes(fill = Group.Size)) + 
  xlab(bquote('Stimulus distance (cm)')) + ylab(bquote('Distance covered (mm)')) + 
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.title=element_blank()) +
  theme(legend.position = "top") 
dc.plot

```

# SCHOOL BEHAVIOR AFTER STIMULATION
Throughout the response, the school’s cohesion and coordination were measured through nearest neighbor distance and alignment, respectively. Nearest neighbor distance (NND) represents the distance between each fish’s center of mass and their most proximal neighbor’s center of mass in the school. Alignment measures the variation in each school’s members orientation with respect to the water's flow (0°). This variation was quantified by using the program Oriana 4 and determining each school members’ angle then calculating the length of mean circular vector (r) of the group, which ranges from 0 (all members are at random angles) to 1 (all angles coordinated exactly). These characteristics were examined at intervals following the stimulus, including 0 ms (representing the school's cohesiveness and coordination immediately prior to the stimulus), 30 ms (representing the typical time for this species to complete stage 1), and 100 ms (the average time for individuals to complete both stages 1 and 2). 

```{r, message = F, warning=FALSE}
#NND
NND1 <- read.csv("/Users/laurennadler/Desktop/MyRData/Monica_groupsize_NND.csv", header=T)
NND1$Groupsize = as.factor(NND1$Groupsize)
NND1$Time = as.factor(NND1$Time)
nnd.lmer = lmer(NND ~  Groupsize*Time + (1|Video/Fish), data=NND1, na.action = "na.omit")
summary(nnd.lmer)
plot(nnd.lmer) #funneled
qqnorm(resid(nnd.lmer)) #curvy
qqline(resid(nnd.lmer))
shapiro.test(resid(nnd.lmer)) #p < 0.001
bartlett.test(resid(nnd.lmer), NND1$Groupsize) #p-value < 0.001

nnd.lmer.bc <- powerTransform(nnd.lmer, family="bcPower") 
summary(nnd.lmer.bc) #both significant

nnd.lmer.bc$roundlam
nnd.bc.lmer = lmer(bcPower(NND, nnd.lmer.bc$roundlam) ~  Groupsize*Time + (1|Video/Fish), data=NND1, na.action = "na.omit")
plot(nnd.bc.lmer) #much better
qqnorm(resid(nnd.bc.lmer)) #much improved
qqline(resid(nnd.bc.lmer)) 
shapiro.test(resid(nnd.bc.lmer)) #p-value = 5.897e-10, but qqplot much improved so continuing with this model
bartlett.test(resid(nnd.bc.lmer), NND1$Groupsize) #p-value = 3.469e-05, but residuals plot much improved so continuing with this model
bartlett.test(resid(nnd.bc.lmer), NND1$Time) #p-value < 2.2e-16, but residuals plot much improved so continuing with this model

nnd.bc.lmer2 = lmer(bcPower(NND, nnd.lmer.bc$roundlam) ~  Groupsize+Time + (1|Video/Fish), data=NND1, na.action = "na.omit")
anova(nnd.bc.lmer, nnd.bc.lmer2) #simpler model lmer2 lower AIC but <2 so keep more complex model
summary(nnd.bc.lmer)
Anova(nnd.bc.lmer, test="F") 
r.squaredGLMM(nnd.bc.lmer)
backtrans.nnd <- update(ref_grid(nnd.bc.lmer), tran = make.tran("boxcox", 0.33)) #where 0.33 is the lambda estimated by powertransform function above
emmeans(backtrans.nnd, pairwise~Groupsize, adjust = c("Tukey"), type = "response") 
emmeans(backtrans.nnd, pairwise~Groupsize*Time, adjust = c("Tukey"), type = "response")

nnd.emm<-emmeans(backtrans.nnd, ~Groupsize*Time, type = "response")
summary(nnd.emm)
nnd.emm<-as.data.frame(nnd.emm)
nnd.emplot = ggplot(nnd.emm, aes(x=Time, y=response, fill=Groupsize)) +
  geom_jitter(aes(x=Time, y = NND, colour=Groupsize), shape = 21, size=1, width=.35, alpha = 0.15, data = NND1) +
  geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL, color = Groupsize), shape=21, size=.2, data = nnd.emm) + 
  geom_line(aes(group = Groupsize, colour = Groupsize), data = nnd.emm) +
  xlab(bquote('Time post-stimulus (ms)')) + ylab(bquote('NND (cm)')) + 
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +  
  theme(axis.text = element_text(size = 12)) + theme(axis.title = element_text(size = 12)) +
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.position = "none") +
  scale_y_log10() 

nnd.emplot + facet_grid(. ~ Groupsize)

#Alignment
align <- read.csv("/Users/laurennadler/Desktop/MyRData/Monica_groupsize_alignment.csv", header=T)
align$Groupsize = as.factor(align$Groupsize)
align$Time = as.factor(align$Time)

alignment.lmer = lmer(Alignment ~ Groupsize*Time + (1|Video), data=align, na.action = "na.omit")
plot(alignment.lmer)
qqnorm(resid(alignment.lmer))
qqline(resid(alignment.lmer))
shapiro.test(resid(alignment.lmer)) #p-value = 0.3076
bartlett.test(resid(alignment.lmer), align$Groupsize) #p-value = 0.6501

alignment.lmer2 = lmer(Alignment ~ Groupsize+Time + (1|Video), data=align, na.action = "na.omit")
anova(alignment.lmer, alignment.lmer2) #similar AIC so keep the more complex model
summary(alignment.lmer)
Anova(alignment.lmer, test="F")
r.squaredGLMM(alignment.lmer)
emmeans(alignment.lmer, pairwise~Time, adjust = c("Tukey"))
emmeans(alignment.lmer, pairwise~Groupsize*Time, adjust = c("Tukey"))

alignment.emm <- ref_grid(alignment.lmer, ~Groupsize*Time, type="response") 
summary(alignment.emm)
alignment.emm<-as.data.frame(alignment.emm)
align.plot = ggplot(alignment.emm, aes(x=Time, y=prediction, fill=Groupsize)) +
  geom_pointrange(aes(ymin = prediction - SE, ymax = prediction + SE, color = Groupsize), na.rm = TRUE, shape=21, size=.7, data = alignment.emm) + 
  geom_line(aes(group = Groupsize, colour = Groupsize), na.rm = TRUE, alignment.emm) +
  geom_jitter(aes(x=Time, y = Alignment, colour=Groupsize), shape = 21, size=1.2, width=.25, alpha = 0.3, data = align) + 
  xlab(bquote('Time post-stimulus (ms)')) + ylab(bquote('Alignment (r)')) +
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +  
  theme(axis.text = element_text(size = 12)) + theme(axis.title = element_text(size = 12)) +
  theme_classic(base_family='Arial', base_size = 28) +
  theme(legend.title=element_blank()) +
  theme(legend.position = "top")
  
align.plot  


```