The influence of spatial context on prey manipulation behaviors in the California moray eel (Gymnothorax mordax)

Research Aims:

Our aim was to assess the impact of spatial context on prey manipulation in the California moray eel (Gymnothorax mordax), a crevice forager with a diverse prey handling repertoire. We first compared the duration engaged in specific manipulation behaviors between spatially “enclosed” and “open” foraging environments (Experiment 1), followed by a comparison of tightly enclosed spaces that were scaled to individual moray diameter (Experiment 2).

Libraries

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(gamlss)
Loading required package: splines
Loading required package: gamlss.data

Attaching package: 'gamlss.data'

The following object is masked from 'package:datasets':

    sleep

Loading required package: gamlss.dist
Loading required package: MASS

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

Loading required package: nlme

Attaching package: 'nlme'

The following object is masked from 'package:dplyr':

    collapse

Loading required package: parallel
 **********   GAMLSS Version 5.4-12  ********** 
For more on GAMLSS look at https://www.gamlss.com/
Type gamlssNews() to see new features/changes/bug fixes.
library(ggpubr)
library(rstatix)

Attaching package: 'rstatix'

The following object is masked from 'package:MASS':

    select

The following object is masked from 'package:stats':

    filter
library(png)
library(lme4)
Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack


Attaching package: 'lme4'

The following object is masked from 'package:gamlss':

    refit

The following object is masked from 'package:nlme':

    lmList
library(ggsignif)
library(glmmTMB)

Attaching package: 'glmmTMB'

The following object is masked from 'package:gamlss':

    refit
library("emmeans")
library(dplyr)
library(MASS)
library(png)
library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:ggpubr':

    get_legend

The following object is masked from 'package:lubridate':

    stamp
library(patchwork)

Attaching package: 'patchwork'

The following object is masked from 'package:cowplot':

    align_plots

The following object is masked from 'package:MASS':

    area
library(magick)
Linking to ImageMagick 6.9.12.93
Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
Disabled features: fftw, ghostscript, x11

EXPERIMENT 1: Open vs Enclosed Crevices

In experiment 1, our objective was to determine whether morays exhibit similar prey manipulation behaviors between open and enclosed foraging environments. We created three experimental feeding treatments; two enclosed (5cm and 10cm in diameter) and one open treatment.

# Naming & reading all the datasets 
P1 <- read_csv(here::here("data/Phase1-SpatialContext.csv")) 
Registered S3 method overwritten by 'bit':
  method   from  
  print.ri gamlss
Rows: 115 Columns: 21
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): Trial info, Treatment, Eel, Knotting_Binomial
dbl (17): Trial_Duration, Handling_Dur, Bites, Ptransport_Dur, Ptransport_Pr...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(P1)
# A tibble: 6 × 21
  `Trial info` Treatment Eel    Trial_Duration Handling_Dur Bites Ptransport_Dur
  <chr>        <chr>     <chr>           <dbl>        <dbl> <dbl>          <dbl>
1 R1 - T1      Open      Bill             19.8        0.855     8           8.84
2 R1 - T1      Open      Johnny           27.1        3.01     12          18.5 
3 R1 - T1      Open      Snick…           77.7       12.3      36          60.6 
4 R1 - T2      Open      Horat…           10.6        0.685     5           5.08
5 R1 - T2      Open      Spagh…           83.2        5.76     24          39   
6 R1 - T2      Open      Eelvis           88         44.5      13           2.8 
# ℹ 14 more variables: Ptransport_Prop <dbl>, Shaking_Dur <dbl>,
#   Shaking_Prop <dbl>, Shaking_Prop_Handling <dbl>, Knotting_Dur <dbl>,
#   Knotting_Prop <dbl>, Knotting_Prop_Handling <dbl>, Knotting_Binomial <chr>,
#   Rotating_Dur <dbl>, Rotating_Prop <dbl>, Rotating_Prop_Handling <dbl>,
#   Ramming_Dur <dbl>, Ramming_Prop <dbl>, Ramming_Prop_Handling <dbl>

P1 data contains: Treatment = ‘Open’ (open spatial treatment), ‘10cm’ (10cm diameter enclosed treatment), ‘5cm’ (5cm diameter enclosed treatment) Eel = name of the individual Trial_Duration = Duration of time in seconds from first bite to completion of swallowing Handling_Dur = Time in seconds spent engaged in at least one of the prey manipulation behaviors Bites = total individual bites Ptransport_Dur = The duration of time in seconds engaged in swallowing (pharyngeal transport)

All behaviors (shaking, ramming, rotating, and knotting) have three columns: Example - Shaking_Dur = the total duration spent shaking Shaking_Prop = Proportion of time spent shaking out of the total Trial_Duration Shaking_Prop_Handling = Proportion of time spent shaking out of the total Handling_Dur

Knotting_Binomial = was knotting observed, yes/no

⌛⌛⌛️ TRIAL COMPONENT DURATIONS ⌛⌛⌛️

This section of code focuses on the total feeding time (Trial_Duration), the time spent handling prey (Handling_Dur), and the time spent swallowing prey (Ptransport_Dur) between treatments (open vs 5cm vs 10cm treatments)

Data Prep - Duration of Trial Components [Feeding Duration, Manipulation Duration, Swallowing Duration]

#Selecting the columns we want:
DurSelP1 <- P1 %>%
  select(Treatment, Eel, Trial_Duration, Ptransport_Dur, Handling_Dur) %>%
  drop_na()

head(DurSelP1) 
# A tibble: 6 × 5
  Treatment Eel        Trial_Duration Ptransport_Dur Handling_Dur
  <chr>     <chr>               <dbl>          <dbl>        <dbl>
1 Open      Bill                 19.8           8.84        0.855
2 Open      Johnny               27.1          18.5         3.01 
3 Open      Snickers             77.7          60.6        12.3  
4 Open      Horatio              10.6           5.08        0.685
5 Open      Spaghettio           83.2          39           5.76 
6 Open      Eelvis               88             2.8        44.5  
#Making treatment a factor:
DurSelP1$Treatment = factor(DurSelP1$Treatment)

## Prepping a dataframe for Graphics - Pivot for long format 
P1_Dur_Long <- DurSelP1 %>%
  pivot_longer(cols = c(Trial_Duration, Ptransport_Dur, Handling_Dur), 
               names_to = "Duration", 
               values_to = "Time")

#Making Duration a factor:
P1_Dur_Long$Duration = factor(P1_Dur_Long$Duration)

head(P1_Dur_Long)
# A tibble: 6 × 4
  Treatment Eel    Duration         Time
  <fct>     <chr>  <fct>           <dbl>
1 Open      Bill   Trial_Duration 19.8  
2 Open      Bill   Ptransport_Dur  8.84 
3 Open      Bill   Handling_Dur    0.855
4 Open      Johnny Trial_Duration 27.1  
5 Open      Johnny Ptransport_Dur 18.5  
6 Open      Johnny Handling_Dur    3.01 
#Calculating Mean & SE For Graphics - This code is used to create future bar graphs comparing the mean time spent engaged in feeding, prey handling, and swallowing between treatments
P1_Dur_SE = P1_Dur_Long %>%
  group_by(Duration, Treatment) %>%
  summarize(mean = mean(Time, na.rm = TRUE),
            sd = sd(Time, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))
`summarise()` has grouped output by 'Duration'. You can override using the
`.groups` argument.
head(P1_Dur_Long)
# A tibble: 6 × 4
  Treatment Eel    Duration         Time
  <fct>     <chr>  <fct>           <dbl>
1 Open      Bill   Trial_Duration 19.8  
2 Open      Bill   Ptransport_Dur  8.84 
3 Open      Bill   Handling_Dur    0.855
4 Open      Johnny Trial_Duration 27.1  
5 Open      Johnny Ptransport_Dur 18.5  
6 Open      Johnny Handling_Dur    3.01 
# --- Getting MEANS & SE for each Individual Duration Category --------------------------

# --------------------------- #
# TOTAL FEEDING DURATION
# --------------------------- #

P1_Feed_SE = DurSelP1 %>%  
  group_by(Treatment) %>% #
  summarize(mean = mean(Trial_Duration, na.rm=TRUE),
            sd = sd(Trial_Duration, na.rm=TRUE),
            n = n(),
            se = sd/sqrt(n))

P1_Feed_SE
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 10cm       30.0  19.1    37  3.13
2 5cm        34.5  23.6    38  3.83
3 Open       38.2  30.4    40  4.81
# --------------------------- #
# HANDLING DURATION
# --------------------------- #

P1_Hand_SE = DurSelP1 %>%  
  group_by(Treatment) %>% #
  summarize(mean = mean(Handling_Dur, na.rm=TRUE),
            sd = sd(Handling_Dur, na.rm=TRUE),
            n = n(),
            se = sd/sqrt(n))

P1_Hand_SE
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 10cm       4.40  6.46    37  1.06
2 5cm        5.29  6.81    38  1.10
3 Open       6.60 14.8     40  2.34
# --------------------------- #
# Swallowing Time
# --------------------------- #

P1_Swallow_SE = DurSelP1 %>%
  group_by(Treatment) %>%
  summarize(mean = mean( Ptransport_Dur, na.rm = TRUE),
            sd = sd(Ptransport_Dur, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))

P1_Swallow_SE 
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 10cm       11.2 11.6     37  1.90
2 5cm        11.6  9.92    38  1.61
3 Open       13.8 12.1     40  1.91

Stats - Duration of Trial Components [Feeding Duration, Manipulation Duration, Swallowing Duration]

This section of code is the statistical analysis comparing each of the Trial Duration componants between the three treatments. For each componant (total time feeding, prey handling, and swallowing), the data is fitted with a GLMM model, with eel identity as a repeated measure…followed by estimated marginal means and pairwise comparisons.

# --------------------------- #
#    TOTAL FEEDING DURATION   #
# --------------------------- #

# Subset the data for "Trial Duration"
TrialDur_Gamma <- DurSelP1 %>%
  select(Eel, Treatment, Trial_Duration) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  # Convert Eel to a factor
         Trial_Duration = as.numeric(Trial_Duration))


# Fit the gamma GLMM with repeated measures
TrialDurationGamma <- glmer(Trial_Duration ~ Treatment + (1 | Eel),
                           data = TrialDur_Gamma,
                           family = Gamma(link = "log"))

# Print model summary
summary(TrialDurationGamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Trial_Duration ~ Treatment + (1 | Eel)
   Data: TrialDur_Gamma

     AIC      BIC   logLik deviance df.resid 
   981.9    995.7   -486.0    971.9      110 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0767 -0.6749 -0.3387  0.2940  3.4881 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 0.09541  0.3089  
 Residual             0.37944  0.6160  
Number of obs: 115, groups:  Eel, 9

Fixed effects:
              Estimate Std. Error t value Pr(>|z|)    
(Intercept)     3.4211     0.1671  20.476   <2e-16 ***
Treatment5cm    0.1390     0.1220   1.140    0.254    
TreatmentOpen   0.1328     0.1245   1.067    0.286    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Trtmn5
Treatmnt5cm -0.366       
TreatmntOpn -0.388  0.498
# Estimated marginal means and pairwise comparisons 
TrialDurationEMMS <- emmeans(TrialDurationGamma, "Treatment")
pairs(TrialDurationEMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm  -0.13902 0.122 Inf  -1.140  0.4897
 10cm - Open -0.13277 0.124 Inf  -1.067  0.5349
 5cm - Open   0.00625 0.124 Inf   0.051  0.9986

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(TrialDurationEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 10cm        3.42 0.167 Inf      3.09      3.75  a    
 Open        3.55 0.165 Inf      3.23      3.88  a    
 5cm         3.56 0.167 Inf      3.23      3.89  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# --------------------------- #
# TOTAL MANIPULATION DURATION #
# --------------------------- #

# Selecting Handling Duration and adding a small constant value to Handling_Dur to handle 0 values
HandDur_Gamma <- DurSelP1 %>%
  select(Eel, Treatment, Handling_Dur) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  
         Handling_Dur = Handling_Dur + 0.001)  # Add a small constant value; rarely, eels would not manipulate prey but would bite and swallow in one movement; to remove 0 values from these isolated instances we added this constant 

# Fit the gamma GLMM with random intercepts
HandDurationGamma<- glmer(Handling_Dur ~ Treatment + (1 | Eel),
                          data = HandDur_Gamma,
                          family = Gamma(link = "log"))

# Print model summary
summary(HandDurationGamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Handling_Dur ~ Treatment + (1 | Eel)
   Data: HandDur_Gamma

     AIC      BIC   logLik deviance df.resid 
   431.9    445.6   -210.9    421.9      110 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.8346 -0.7790 -0.3567  0.3871  3.6862 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 1.734    1.317   
 Residual             1.435    1.198   
Number of obs: 115, groups:  Eel, 9

Fixed effects:
              Estimate Std. Error t value Pr(>|z|)   
(Intercept)     1.2392     0.4450   2.785  0.00535 **
Treatment5cm    0.1383     0.3741   0.370  0.71159   
TreatmentOpen  -0.1137     0.3940  -0.289  0.77295   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Trtmn5
Treatmnt5cm -0.416       
TreatmntOpn -0.436  0.467
# Estimated marginal means and pairwise comparisons 
HandDurationEMMS <- emmeans(HandDurationGamma, "Treatment")
pairs(HandDurationEMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm    -0.138 0.374 Inf  -0.370  0.9274
 10cm - Open    0.114 0.394 Inf   0.289  0.9551
 5cm - Open     0.252 0.397 Inf   0.635  0.8010

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(HandDurationEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Open        1.13 0.448 Inf     0.248      2.00  a    
 10cm        1.24 0.445 Inf     0.367      2.11  a    
 5cm         1.38 0.447 Inf     0.502      2.25  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# --------------------------- #
#     Pharyngeal Transport    #
# --------------------------- #

# Selecting Pharyngeal Transport duration and adding a small constant value to  Ptransport_Dur to handle 0 values
PT_Gamma <- DurSelP1 %>%
  select(Eel, Treatment, Ptransport_Dur) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  
         Ptransport_Dur = Ptransport_Dur+ 0.001)  # Add a small constant value

# Fit the gamma GLMM with random intercepts
PTGamma<- glmer(Ptransport_Dur ~ Treatment + (1 | Eel),
                          data = PT_Gamma,
                          family = Gamma(link = "log"))
boundary (singular) fit: see help('isSingular')
# Print model summary
summary(PTGamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Ptransport_Dur ~ Treatment + (1 | Eel)
   Data: PT_Gamma

     AIC      BIC   logLik deviance df.resid 
   817.8    831.5   -403.9    807.8      110 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.1001 -0.6842 -0.3837  0.2845  3.7468 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 Eel      (Intercept) 1.345e-14 1.160e-07
 Residual             8.261e-01 9.089e-01
Number of obs: 115, groups:  Eel, 9

Fixed effects:
              Estimate Std. Error t value Pr(>|z|)    
(Intercept)    2.42033    0.17399  13.910   <2e-16 ***
Treatment5cm   0.03414    0.24444   0.140    0.889    
TreatmentOpen  0.20113    0.24141   0.833    0.405    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Trtmn5
Treatmnt5cm -0.712       
TreatmntOpn -0.721  0.513
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
# Estimated marginal means and pairwise comparisons 
PT_EMMS <- emmeans(PTGamma, "Treatment")
pairs(PT_EMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm   -0.0341 0.244 Inf  -0.140  0.9893
 10cm - Open  -0.2011 0.241 Inf  -0.833  0.6824
 5cm - Open   -0.1670 0.240 Inf  -0.697  0.7655

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(PT_EMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 10cm        2.42 0.174 Inf      2.08      2.76  a    
 5cm         2.45 0.172 Inf      2.12      2.79  a    
 Open        2.62 0.167 Inf      2.29      2.95  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Graphics - Durations of Trial Components

The next code makes the bar graphs for the trial duration components, comparing the durations feeding, handling prey, and swallowing between the three treatments (open, 10cm, 5cm)

#Naming the Future Facets: 
Facet_Names_P1_Dur <- list(
  'Trial_Duration' = "Total Feeding Duration",
  'Handling_Dur' = "Prey Manipulation Duration",
  'Ptransport_Dur' = "Swallowing Duration")

Facet_labeller_P1_Dur <- function(variable,value){
  return(Facet_Names_P1_Dur[value])
}

#The facets of the graph will be the different duration componants: Feeding, Handling, Swallowing

#Making a plot
P1_Dur_Graph = P1_Dur_SE %>% 
  mutate(Treatment = fct_relevel(Treatment, "5cm", "10cm", "Open")) %>% 
  mutate(Duration = fct_relevel(Duration, "Trial_Duration","Handling_Dur", "Ptransport_Dur")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           position = position_dodge()) +
  # error bar 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .25) + 
  scale_fill_manual(values = c("darkslategray4","darkslategray3", "azure2")) +
  #scale_x_discrete(labels=c("Small","Medium", "Large"))+
  labs(y = "Mean Duration (Sec)", x = "") +
  # facet wrap + the label code from above for our new label names
  facet_wrap(. ~ Duration, labeller = Facet_labeller_P1_Dur, scales = "free") +
  theme_classic() +
  theme (axis.text.x = element_text(size = 15),
        axis.text.y =  element_text(size = 15),
        axis.title.y = element_text(size = 20),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        plot.title = element_text(hjust =.5, size = 10, colour = "black"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face="bold"))
Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
arguments are now deprecated.
ℹ See labellers documentation.
P1_Dur_Graph

# Saves to reports folder within the project
ggsave(filename = "P1-TrialComponents.png", plot = P1_Dur_Graph, width = 10, height = 5, dpi = 500)

⏰⏰⏰ BEHAVIOR DURATION ⏰⏰⏰

This next chapter of code compares the time spent engaged in different prey manipulation behaviors between treatments. Because there was no significant difference in total manipulation time (see previous) between treatments, we used the total duration engaged in each behavior (Example Column: Shaking_Duration), as opposed to proportional data (out of total feeding or handling time)

Data Prep - Duration of Behaviors

head(P1)
# A tibble: 6 × 21
  `Trial info` Treatment Eel    Trial_Duration Handling_Dur Bites Ptransport_Dur
  <chr>        <chr>     <chr>           <dbl>        <dbl> <dbl>          <dbl>
1 R1 - T1      Open      Bill             19.8        0.855     8           8.84
2 R1 - T1      Open      Johnny           27.1        3.01     12          18.5 
3 R1 - T1      Open      Snick…           77.7       12.3      36          60.6 
4 R1 - T2      Open      Horat…           10.6        0.685     5           5.08
5 R1 - T2      Open      Spagh…           83.2        5.76     24          39   
6 R1 - T2      Open      Eelvis           88         44.5      13           2.8 
# ℹ 14 more variables: Ptransport_Prop <dbl>, Shaking_Dur <dbl>,
#   Shaking_Prop <dbl>, Shaking_Prop_Handling <dbl>, Knotting_Dur <dbl>,
#   Knotting_Prop <dbl>, Knotting_Prop_Handling <dbl>, Knotting_Binomial <chr>,
#   Rotating_Dur <dbl>, Rotating_Prop <dbl>, Rotating_Prop_Handling <dbl>,
#   Ramming_Dur <dbl>, Ramming_Prop <dbl>, Ramming_Prop_Handling <dbl>
#Selecting the columns we want:
SelBehaviorDurationsP1 <- P1 %>%
  select(Treatment, Eel, Shaking_Dur, Ramming_Dur, Knotting_Dur, Rotating_Dur) %>%
  drop_na()

head(SelBehaviorDurationsP1) 
# A tibble: 6 × 6
  Treatment Eel        Shaking_Dur Ramming_Dur Knotting_Dur Rotating_Dur
  <chr>     <chr>            <dbl>       <dbl>        <dbl>        <dbl>
1 Open      Bill             0.855        0            0           0    
2 Open      Johnny           0            0            0           3.01 
3 Open      Snickers         0.97         1.75         6.5         3.06 
4 Open      Horatio          0.685        0            0           0    
5 Open      Spaghettio       0            0            4.46        1.3  
6 Open      Eelvis           0            7.7         36.3         0.518
# Changing to long format 
P1BehaviorDur_Long <- SelBehaviorDurationsP1 %>%
  pivot_longer(cols = c(Shaking_Dur, Ramming_Dur, Knotting_Dur, Rotating_Dur), 
               names_to = "Behavior", 
               values_to = "Duration")


#Calculating Mean & SE For Graphics
P1BehaviorDurSE <- P1BehaviorDur_Long %>%
  group_by(Treatment, Behavior) %>%
  summarize(mean = mean(Duration, na.rm = TRUE),
            sd = sd(Duration, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))
`summarise()` has grouped output by 'Treatment'. You can override using the
`.groups` argument.
P1BehaviorDurSE
# A tibble: 12 × 6
# Groups:   Treatment [3]
   Treatment Behavior       mean     sd     n     se
   <chr>     <chr>         <dbl>  <dbl> <int>  <dbl>
 1 10cm      Knotting_Dur 1.76    2.92     37 0.480 
 2 10cm      Ramming_Dur  0.0954  0.390    37 0.0642
 3 10cm      Rotating_Dur 1.92    4.37     37 0.718 
 4 10cm      Shaking_Dur  0.622   0.759    37 0.125 
 5 5cm       Knotting_Dur 1.47    3.91     38 0.635 
 6 5cm       Ramming_Dur  0.228   1.17     38 0.190 
 7 5cm       Rotating_Dur 2.92    4.71     38 0.764 
 8 5cm       Shaking_Dur  0.676   0.845    38 0.137 
 9 Open      Knotting_Dur 4.62   12.2      40 1.93  
10 Open      Ramming_Dur  0.877   2.77     40 0.438 
11 Open      Rotating_Dur 0.692   0.940    40 0.149 
12 Open      Shaking_Dur  0.403   0.579    40 0.0916
#Making treatment a factor:
P1BehaviorDurSE$Treatment = factor(P1BehaviorDurSE$Treatment)

Stats - Duration of Behaviors

The next code is the statistical analysis, comparing the time spent engaged in each behavior (ramming, rotating, shaking, knotting) between the treatments (open, 10cm, 5cm). To compare average duration spent engaged in each behavior, we fit the data for each behavior with a GLMM using a Tweedie distribution. The data is not normal and shows overdisperson making other distributions less suitable. The Tweedie distribution encompasses a wide range of distribution shapes and is also well suited for modeling continuous positive data inclusive of 0 values.

Each behavior has an emoji header to quickly find each chunk.

# 😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫
#  ROTATING STATS
# 😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫


# Pulling out the ROTATING duration data
RotatingTweedie <- P1BehaviorDur_Long %>% 
  filter(Behavior == "Rotating_Dur")

# Fit the model using glmmTMB
RotatingTweedie_Mod  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                    data = RotatingTweedie,
                                    family = tweedie)

# Print the model summary
summary(RotatingTweedie_Mod)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: RotatingTweedie

     AIC      BIC   logLik deviance df.resid 
   422.9    439.3   -205.4    410.9      109 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.5377   0.7333  
Number of obs: 115, groups:  Eel, 9

Dispersion parameter for tweedie family (): 2.88 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)     0.4285     0.3564   1.202   0.2292  
Treatment5cm    0.4476     0.3196   1.400   0.1614  
TreatmentOpen  -0.9017     0.3894  -2.316   0.0206 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
RotatingDurEMMS <- emmeans(RotatingTweedie_Mod, "Treatment")
pairs(RotatingDurEMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm    -0.448 0.320 Inf  -1.400  0.3407
 10cm - Open    0.902 0.389 Inf   2.316  0.0537
 5cm - Open     1.349 0.371 Inf   3.632  0.0008

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(RotatingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Open      -0.473 0.392 Inf    -1.241     0.295  a    
 10cm       0.429 0.356 Inf    -0.270     1.127  ab   
 5cm        0.876 0.336 Inf     0.219     1.534   b   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# 🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏
#   RAMMING STATS
# 🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏


# Pulling out the RAMMING duration data
RammingTweedie <- P1BehaviorDur_Long %>% 
  filter(Behavior == "Ramming_Dur")

# Fit the model using glmmTMB
RammingTweedie_Mod  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                  data = RammingTweedie,
                                  family = tweedie)

# Print the model summary
summary(RammingTweedie_Mod)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: RammingTweedie

     AIC      BIC   logLik deviance df.resid 
   139.1    155.6    -63.6    127.1      109 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 1.591    1.261   
Number of obs: 115, groups:  Eel, 9

Dispersion parameter for tweedie family ():  7.3 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)   
(Intercept)    -2.7457     1.0629  -2.583  0.00978 **
Treatment5cm    0.6843     1.1563   0.592  0.55397   
TreatmentOpen   1.8185     1.0700   1.700  0.08922 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
RammingDurEMMS <- emmeans(RammingTweedie_Mod, "Treatment")
pairs(RammingDurEMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm    -0.684 1.156 Inf  -0.592  0.8245
 10cm - Open   -1.819 1.070 Inf  -1.700  0.2053
 5cm - Open    -1.134 0.951 Inf  -1.193  0.4571

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(RammingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 10cm      -2.746 1.063 Inf     -4.83    -0.663  a    
 5cm       -2.061 0.985 Inf     -3.99    -0.131  a    
 Open      -0.927 0.878 Inf     -2.65     0.793  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# 🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝
#  SHAKING STATS 
# 🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝

# Pulling out the SHAKING DURATION Data
ShakingTweedie <- P1BehaviorDur_Long %>% 
  filter(Behavior == "Shaking_Dur")


# Fit the model using glmmTMB
ShakingTweedie_Mod  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                 data = ShakingTweedie,
                                 family = tweedie)

# Print the model summary
summary(ShakingTweedie_Mod)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: ShakingTweedie

     AIC      BIC   logLik deviance df.resid 
   246.5    263.0   -117.2    234.5      109 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.9616   0.9806  
Number of obs: 115, groups:  Eel, 9

Dispersion parameter for tweedie family (): 0.744 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)    -0.7446     0.3813  -1.953   0.0508 .
Treatment5cm    0.1049     0.2446   0.429   0.6680  
TreatmentOpen  -0.5539     0.2712  -2.042   0.0411 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
ShakingDurEMMS <- emmeans(ShakingTweedie_Mod, "Treatment")
pairs(ShakingDurEMMS)
 contrast    estimate    SE  df z.ratio p.value
 10cm - 5cm    -0.105 0.245 Inf  -0.429  0.9036
 10cm - Open    0.554 0.271 Inf   2.042  0.1022
 5cm - Open     0.659 0.267 Inf   2.464  0.0366

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(ShakingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Open      -1.298 0.397 Inf     -2.08  -0.51943  a    
 10cm      -0.745 0.381 Inf     -1.49   0.00269  ab   
 5cm       -0.640 0.378 Inf     -1.38   0.10073   b   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
#🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢
# KNOTTING STATS 
#🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢🪢

#Making treatment a factor:
P1BehaviorDur_Long$Treatment = factor(P1BehaviorDur_Long$Treatment)


# Pulling out the SHAKING DURATION Data
KnottingTweedie <- P1BehaviorDur_Long %>% 
  filter(Behavior == "Knotting_Dur") %>% 
   mutate(Treatment = relevel(Treatment, ref = "Open"))


# Fit the model using glmmTMB
KnottingTweedie_Mod  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                 data = KnottingTweedie,
                                 family = tweedie)

# Print the model summary
summary(KnottingTweedie_Mod)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: KnottingTweedie

     AIC      BIC   logLik deviance df.resid 
   341.4    357.9   -164.7    329.4      109 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 2.505    1.583   
Number of obs: 115, groups:  Eel, 9

Dispersion parameter for tweedie family ():  5.4 

Conditional model:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)     0.4291     0.6533   0.657   0.5113  
Treatment10cm  -0.5803     0.4129  -1.406   0.1599  
Treatment5cm   -0.7944     0.4450  -1.785   0.0743 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
KnottingDurEMMS <- emmeans(KnottingTweedie_Mod, "Treatment")
pairs(KnottingDurEMMS)
 contrast    estimate    SE  df z.ratio p.value
 Open - 10cm    0.580 0.413 Inf   1.405  0.3380
 Open - 5cm     0.794 0.445 Inf   1.785  0.1746
 10cm - 5cm     0.214 0.488 Inf   0.438  0.8996

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(KnottingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 5cm       -0.365 0.695 Inf    -1.728     0.998  a    
 10cm      -0.151 0.671 Inf    -1.467     1.165  a    
 Open       0.429 0.653 Inf    -0.851     1.710  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Graphics - Duration of Behaviors

# This graph shows the average DURATION spent in each behavior.
# 0's have been left in the dataset 

#Naming the Future Facets:
Facet_Names_P1_BehaviorDur <- list(
  'Knotting_Dur' = "KNOTTING",
  'Ramming_Dur' = "RAMMING",
  'Rotating_Dur' = "ROTATING",
  'Shaking_Dur' = "SHAKING")

Facet_labeller_BehaviorDur<- function(variable,value){
  return(Facet_Names_P1_BehaviorDur[value])}

#Wrangling some pictures for the future
KnotPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Knotting.png") 
ShakePic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Shaking.png")
RotatePic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Rotating.png")
RamPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Ramming.png")

# Creating a data frame with pairwise significance: 
p1_handdur_signif <- tribble(
  ~Treatment, ~Behavior,     ~Signif,
  "Open",   "Knotting_Dur", "a",    
  "10cm",   "Knotting_Dur", "a",   
  "5cm",    "Knotting_Dur", "a",    
  "Open",    "Ramming_Dur", "a",    
  "10cm",   "Ramming_Dur", "a",   
  "5cm",    "Ramming_Dur", "a",   
  "Open",    "Rotating_Dur", "b",    
  "10cm",   "Rotating_Dur", "ab",   
  "5cm",    "Rotating_Dur", "a",    
  "Open",    "Shaking_Dur", "b",    
  "10cm",   "Shaking_Dur", "ab",   
  "5cm",    "Shaking_Dur", "a" 
) %>% 
  left_join(transmute(P1BehaviorDurSE, Treatment, Behavior, mean, se),
            by = c("Treatment", "Behavior")) %>% 
  mutate(Treatment = fct_relevel(Treatment, "5cm", "10cm", "Open"),
         Behavior = factor(Behavior),
         y = mean + se + 0.5)


# The graph:
P1_BehaviorDur = P1BehaviorDurSE %>% 
  mutate(Treatment = fct_relevel(Treatment, "5cm", "10cm", "Open")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           position = position_dodge()) +
  #Adding error bars: 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .5) + 
  #Adding Significance Labels
  geom_text(aes(y = y, label = Signif), p1_handdur_signif, size =5) +
  #All other aesthetics:
  scale_fill_manual(values = c("darkslategray4", "darkslategray3", "azure2")) +
  scale_x_discrete(labels=c("5cm","10cm", "Open"))+
  labs(y = "Mean Duration in Behavior (Sec)", x = "Treatment") +
  facet_wrap(. ~ Behavior, labeller = Facet_labeller_BehaviorDur) +
  theme_classic() +
  theme (axis.text.x = element_text(size = 15),
         axis.text.y =  element_text(size = 15),
         axis.title.x = element_text(size = 20),
         axis.title.y = element_text(size = 20),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        plot.title = element_text(hjust =.5, size = 40, colour = "black"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face = "bold"))
Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
arguments are now deprecated.
ℹ See labellers documentation.
P1_BehaviorDur

P1_BehaviorDur_Pics <- ggdraw(P1_BehaviorDur) +
draw_image(KnotPic, x=0.35 , y=0.89, hjust = 1, vjust =1, height = .17, width
=.17) +
draw_image(RamPic, x=0.85 , y=0.89, hjust = 1, vjust =1, height = .2, width
=.2 ) +
draw_image(RotatePic, x=0.43 , y=0.5, hjust = 1, vjust =1, height = .24, width
=.24 ) +
draw_image(ShakePic, x=0.86 , y=0.46, hjust = 1, vjust =1, height = .2, width
=.2 ) 

P1_BehaviorDur_Pics 

ggsave(filename = "P1-BehaviorDuration.png", plot = P1_BehaviorDur_Pics, width = 10, height = 5, dpi = 300)

————————————

Experiment 1 - FINAL FIGURE Graphic

———————————–

The code below combines the Trial Component Figure with the Behavior Duration Figure for Experiment 1

P1_FinalFigure <- ggarrange(P1_Dur_Graph, P1_BehaviorDur_Pics,
                    labels = c("A", "B"),
                    ncol = 1, nrow = 2,
                    heights = c(0.25, 0.5))
P1_FinalFigure

ggsave(filename = "P1_Final.png", plot = P1_FinalFigure, width = 11, height = 10, dpi = 300)

EXPERIMENT 2: Enclosed Spatial Trials

In Experiment 2, our aim was to document the spatial threshold relative to eel size that might limit more three-dimensional movement (i.e. knotting). Treatments (small, medium, and large) were scaled to individual diameters of the moray. The large treatment was 2 x the eel diameter, the medium treatment 1.5 x the eel’s diameter,and the small treatment we scaled to be the exact diameter of the eel plus 10% of the eel’s diameter.

We included two new behaviors in our analysis: tail and body anchoring

#Reading in the dataset
P2 <- read_csv(here::here("data/Phase2-SpatialContext2.csv"))
Rows: 58 Columns: 37
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): Eel, Treatment, Knotting_Binomial
dbl (34): Treatment_Round, Trial_Duration, Handling_Dur, Handling-PT, PT_Dur...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(P2) 
# A tibble: 6 × 37
  Eel        Treatment Treatment_Round Trial_Duration Handling_Dur `Handling-PT`
  <chr>      <chr>               <dbl>          <dbl>        <dbl>         <dbl>
1 Eelvis     Small                   1           50.4        16.1           49.3
2 Snickers   Small                   1           32.2         8.05          29.6
3 Bill       Small                   1           27.4        10.3           19.4
4 Hamlet     Small                   1           38.2        15.0           24.8
5 Horatio    Small                   1           30.4        40.7           17.9
6 Spaghettio Small                   1           45.9        51.7           42.9
# ℹ 31 more variables: PT_Duration <dbl>, PT_Prop <dbl>,
#   Knotting_Duration <dbl>, Knotting_Prop <dbl>, Knotting_Prop_Handling <dbl>,
#   Knotting_Binomial <chr>, Rotating_Duration <dbl>, Rotating_Prop <dbl>,
#   Rotating_Prop_Handling <dbl>, Ramming_Duration <dbl>, Ramming_Prop <dbl>,
#   Ramming_Prop_Handling <dbl>, Shaking_Duration <dbl>, Shaking_Prop <dbl>,
#   Shaking_Prop_Handling <dbl>, Crinkle_Duration <dbl>, Crinkle_Prop <dbl>,
#   Crinkle_Prop_Handling <dbl>, Tug_Duration <dbl>, Tug_Prop <dbl>, …

P2 data contains: Treatment = ‘small’, ‘medium’, ‘large’ Eel = name of the individual Treatment_Round = individual’s exposure to the treatment (1st time, 2nd time, etc.) Trial_Duration = Duration of time in seconds from first bite to completion of swallowing Handling_Duration = Time in seconds spent engaged in at least one of the prey manipulation behaviors PT_Duration = The duration of time in seconds engaged in swallowing (pharyngeal transport)

All behaviors (shaking, ramming, rotating, knotting, crinkle, anchor, tug) have three columns: Example - Shaking_Dur = the total duration spent shaking Shaking_Prop = Proportion of time spent shaking out of the total Trial_Duration Shaking_Prop_Handling = Proportion of time spent shaking out of the total Handling_Dur

  • Tug was not included in analysis, as it could not be fully seperated from other behaviors
  • ‘crinkle’ in the data set is code for ‘body anchoring’ in the manuscript; where an eel forms a ‘w’ shape, pressing on the tube walls.
  • ‘anchor’ in the data set refers specifically to tail anchoring

Knotting_Binomial = was knotting observed, yes/no

⌛️⌛️⌛️ TRIAL COMPONENT DURATION ⌛⌛⌛️

This section of code focuses on the total feeding time (Trial_Duration), the time spent handling prey (Handling_Duration), and the time spent swallowing prey (PT_Duration) between treatments (small, medium, large)

Data Prep - Duration of Trial Components

#Selecting the columns we want:
DurSelP2 <- P2 %>%
  select(Treatment, Eel, Trial_Duration, PT_Duration, Handling_Dur) %>%
  drop_na()

head(DurSelP2) 
# A tibble: 6 × 5
  Treatment Eel        Trial_Duration PT_Duration Handling_Dur
  <chr>     <chr>               <dbl>       <dbl>        <dbl>
1 Small     Eelvis               50.4        1.17        16.1 
2 Small     Snickers             32.2        2.62         8.05
3 Small     Bill                 27.4        7.99        10.3 
4 Small     Hamlet               38.2       13.3         15.0 
5 Small     Horatio              30.4       12.5         40.7 
6 Small     Spaghettio           45.9        3.01        51.7 
#Making treatment a factor:
DurSelP2$Treatment = factor(DurSelP2$Treatment)

## Prepping a dataframe for Graphics 
P2_Dur_Long <- DurSelP2 %>%
  pivot_longer(cols = c(Trial_Duration, PT_Duration, Handling_Dur), 
               names_to = "Duration", 
               values_to = "Time")

#Making Duration a factor:
P2_Dur_Long$Duration = factor(P2_Dur_Long$Duration)

head(P2_Dur_Long)
# A tibble: 6 × 4
  Treatment Eel      Duration        Time
  <fct>     <chr>    <fct>          <dbl>
1 Small     Eelvis   Trial_Duration 50.4 
2 Small     Eelvis   PT_Duration     1.17
3 Small     Eelvis   Handling_Dur   16.1 
4 Small     Snickers Trial_Duration 32.2 
5 Small     Snickers PT_Duration     2.62
6 Small     Snickers Handling_Dur    8.05
## Calculating Mean & SE For Graphics - This code is used to create future bar graphs comparing the mean time spent engaged in feeding, prey handling, and swallowing between treatments
P2_Dur_SE = P2_Dur_Long %>%
  group_by(Duration, Treatment) %>%
  summarize(mean = mean(Time, na.rm = TRUE),
            sd = sd(Time, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))
`summarise()` has grouped output by 'Duration'. You can override using the
`.groups` argument.
head(P2_Dur_SE)
# A tibble: 6 × 6
# Groups:   Duration [2]
  Duration     Treatment  mean    sd     n    se
  <fct>        <fct>     <dbl> <dbl> <int> <dbl>
1 Handling_Dur Large     16.6   9.66    21  2.11
2 Handling_Dur Medium    21.4  17.0     20  3.81
3 Handling_Dur Small     23.8  15.1     17  3.65
4 PT_Duration  Large      8.41  6.14    21  1.34
5 PT_Duration  Medium     8.64  6.01    20  1.34
6 PT_Duration  Small      8.72  7.29    17  1.77
# --- Getting MEANS & SE for each Individual Duration Category --------------------------

# --------------------------- #
# TOTAL FEEDING DURATION
# --------------------------- #

P2_Feed_SE=DurSelP2 %>%  
  group_by(Treatment) %>% #
  summarize(mean=mean(Trial_Duration, na.rm=TRUE),
            sd=sd(Trial_Duration, na.rm=TRUE),
            n=n(),
            se=sd/sqrt(n))

P2_Feed_SE
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 Large      47.3  35.3    21  7.69
2 Medium     62.2  68.4    20 15.3 
3 Small      39.4  17.0    17  4.13
# --------------------------- #
# HANDLING DURATION
# --------------------------- #

P2_Hand_SE=DurSelP2%>%  
  group_by(Treatment) %>% #
  summarize(mean=mean(Handling_Dur, na.rm=TRUE),
            sd=sd(Handling_Dur, na.rm=TRUE),
            n=n(),
            se=sd/sqrt(n))

P2_Hand_SE
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 Large      16.6  9.66    21  2.11
2 Medium     21.4 17.0     20  3.81
3 Small      23.8 15.1     17  3.65
# --------------------------- #
# SWALLOWING DURATION
# --------------------------- #

P2_Swallow_SE = DurSelP2 %>%
  group_by(Treatment) %>%
  summarize(mean = mean( PT_Duration, na.rm = TRUE),
            sd = sd(PT_Duration, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))

P2_Swallow_SE 
# A tibble: 3 × 5
  Treatment  mean    sd     n    se
  <fct>     <dbl> <dbl> <int> <dbl>
1 Large      8.41  6.14    21  1.34
2 Medium     8.64  6.01    20  1.34
3 Small      8.72  7.29    17  1.77

Stats - Duration of Trial Components

This section of code is the statistical analysis comparing each of the Trial Duration componants between the three treatments. For each componant (total time feeding, prey handling, and swallowing), the data is fitted with a GLMM model, with eel identity as a repeated measure…followed by estimated marginal means and pairwise comparisons.

# Subset the data for "Trial Duration"
TrialDur2_Gamma <- DurSelP2 %>%
  select(Eel, Treatment, Trial_Duration) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  # Convert Eel to a factor
         Trial_Duration = as.numeric(Trial_Duration))


# Fit the gamma GLMM with repeated measures, adjusting optimization settings
TrialDuration2Gamma <- glmer(Trial_Duration ~ Treatment + (1 | Eel),
                           data = TrialDur2_Gamma,
                           family = Gamma(link = "log"))
 
# Print model summary
summary(TrialDuration2Gamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Trial_Duration ~ Treatment + (1 | Eel)
   Data: TrialDur2_Gamma

     AIC      BIC   logLik deviance df.resid 
   546.6    556.9   -268.3    536.6       53 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.0527 -0.7165 -0.2225  0.4843  3.2157 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 0.1528   0.3908  
 Residual             0.3934   0.6272  
Number of obs: 58, groups:  Eel, 8

Fixed effects:
                Estimate Std. Error t value Pr(>|z|)    
(Intercept)      3.80443    0.22848  16.651   <2e-16 ***
TreatmentMedium  0.03748    0.17664   0.212    0.832    
TreatmentSmall  -0.03303    0.17636  -0.187    0.851    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) TrtmnM
TreatmntMdm -0.331       
TretmntSmll -0.325  0.419
TrialDuration2EMMS <- emmeans(TrialDuration2Gamma, "Treatment")
pairs(TrialDuration2EMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium  -0.0375 0.177 Inf  -0.212  0.9755
 Large - Small    0.0330 0.176 Inf   0.187  0.9808
 Medium - Small   0.0705 0.190 Inf   0.371  0.9271

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(TrialDuration2EMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Small       3.77 0.239 Inf      3.30      4.24  a    
 Large       3.80 0.228 Inf      3.36      4.25  a    
 Medium      3.84 0.238 Inf      3.38      4.31  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# Handling Duration  #

# Add a small constant value to Handling_Dur to handle 0 values
HandDur2_Gamma <- DurSelP2 %>%
  select(Eel, Treatment, Handling_Dur) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  
         Handling_Dur = Handling_Dur + 0.001)  # Add a small constant value

# Fit the gamma GLMM with random intercepts
HandDuration2Gamma<- glmer(Handling_Dur ~ Treatment + (1 | Eel),
                          data = HandDur2_Gamma,
                          family = Gamma(link = "log"))

# Print model summary
summary(HandDuration2Gamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Handling_Dur ~ Treatment + (1 | Eel)
   Data: HandDur2_Gamma

     AIC      BIC   logLik deviance df.resid 
   462.6    472.9   -226.3    452.6       53 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.4526 -0.7695 -0.3683  0.6058  3.0399 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 0.002535 0.05035 
 Residual             0.435974 0.66028 
Number of obs: 58, groups:  Eel, 8

Fixed effects:
                Estimate Std. Error t value Pr(>|z|)    
(Intercept)       2.8065     0.1557  18.028   <2e-16 ***
TreatmentMedium   0.2499     0.2220   1.125    0.260    
TreatmentSmall    0.3640     0.2295   1.586    0.113    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) TrtmnM
TreatmntMdm -0.655       
TretmntSmll -0.662  0.442
HandDuration2EMMS <- emmeans(HandDuration2Gamma, "Treatment")
pairs(HandDuration2EMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium   -0.250 0.222 Inf  -1.125  0.4984
 Large - Small    -0.364 0.229 Inf  -1.586  0.2516
 Medium - Small   -0.114 0.238 Inf  -0.479  0.8814

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(HandDuration2EMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Large       2.81 0.156 Inf      2.50      3.11  a    
 Medium      3.06 0.168 Inf      2.73      3.39  a    
 Small       3.17 0.172 Inf      2.83      3.51  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# Pharyngeal Transport #

PT2_Gamma <- DurSelP2 %>%
  select(Eel, Treatment, PT_Duration) %>%
  mutate(Treatment = factor(Treatment),
         Eel = factor(Eel),  
         PT_Duration = PT_Duration + 0.001)  # Add a small constant value

# Fit the gamma GLMM with random intercepts
PT2Gamma<- glmer(PT_Duration ~ Treatment + (1 | Eel),
                          data = PT2_Gamma,
                          family = Gamma(link = "log"))
boundary (singular) fit: see help('isSingular')
# Print model summary
summary(PT2Gamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: PT_Duration ~ Treatment + (1 | Eel)
   Data: PT2_Gamma

     AIC      BIC   logLik deviance df.resid 
   377.8    388.1   -183.9    367.8       53 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3674 -0.8195 -0.1740  0.5566  3.5107 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 0.0000   0.0000  
 Residual             0.5347   0.7312  
Number of obs: 58, groups:  Eel, 8

Fixed effects:
                Estimate Std. Error t value Pr(>|z|)    
(Intercept)      2.12949    0.22764   9.355   <2e-16 ***
TreatmentMedium  0.02703    0.32593   0.083    0.934    
TreatmentSmall   0.03678    0.34034   0.108    0.914    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) TrtmnM
TreatmntMdm -0.698       
TretmntSmll -0.669  0.467
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
PT2_EMMS <- emmeans(PT2Gamma, "Treatment")
pairs(PT2_EMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium -0.02703 0.326 Inf  -0.083  0.9962
 Large - Small  -0.03678 0.340 Inf  -0.108  0.9936
 Medium - Small -0.00976 0.344 Inf  -0.028  0.9996

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(PT2_EMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Large       2.13 0.228 Inf      1.68      2.58  a    
 Medium      2.16 0.233 Inf      1.70      2.61  a    
 Small       2.17 0.253 Inf      1.67      2.66  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Graphics - Duration of Trial Components

The next code makes the bar graphs for the trial duration components, comparing the durations feeding, handling prey, and swallowing between the three treatments (small, medium, large).

#GRAPH FACETS = DURATIONS 

#Naming the Future Facets
Facet_Names_P2_Dur <- list(
  'Trial_Duration' = "Total Feeding Duration",
  'Handling_Dur' = "Prey Manipulation Duration",
  'PT_Duration' = "Swallowing Duration")

Facet_labeller_P2_Dur <- function(variable,value){
  return(Facet_Names_P2_Dur[value])
}


#Making a plot

P2_Dur_Graph = P2_Dur_SE %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large")) %>% 
  mutate(Duration = fct_relevel(Duration, "Trial_Duration","Handling_Dur", "PT_Duration")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           position = position_dodge()) +
  # error bar 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .25) + 
  scale_fill_manual(values = c("darkslategray4","darkslategray3", "azure2")) +
  labs(y = "Mean Duration (Sec)", x = "") +
  # facet wrap + the label code from above for our new label names
  facet_wrap(. ~ Duration, labeller = Facet_labeller_P2_Dur, scales = "free") +
  theme_classic() +
theme (axis.text.x = element_text(size = 15),
        axis.text.y =  element_text(size = 15),
        axis.title.y = element_text(size = 20),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face = "bold"))
Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
arguments are now deprecated.
ℹ See labellers documentation.
P2_Dur_Graph

# Saves to reports folder within the project
ggsave(filename = "P2-Trial Components.png", plot = P2_Dur_Graph, width = 11, height = 5, dpi = 300)

⏰ ⏰ ⏰ BEHAVIOR DURATION ⏰⏰⏰

This next chapter of code compares the time spent engaged in different prey manipulation behaviors between treatments. Because there was no significant difference in total manipulation time (see previous) between treatments, we used the total duration engaged in each behavior (Example Column: ‘Shaking_Duration’), as opposed to proportional data (out of total feeding or handling time)

Data Prep - Duration of Behaviors

head(P2)
# A tibble: 6 × 37
  Eel        Treatment Treatment_Round Trial_Duration Handling_Dur `Handling-PT`
  <chr>      <chr>               <dbl>          <dbl>        <dbl>         <dbl>
1 Eelvis     Small                   1           50.4        16.1           49.3
2 Snickers   Small                   1           32.2         8.05          29.6
3 Bill       Small                   1           27.4        10.3           19.4
4 Hamlet     Small                   1           38.2        15.0           24.8
5 Horatio    Small                   1           30.4        40.7           17.9
6 Spaghettio Small                   1           45.9        51.7           42.9
# ℹ 31 more variables: PT_Duration <dbl>, PT_Prop <dbl>,
#   Knotting_Duration <dbl>, Knotting_Prop <dbl>, Knotting_Prop_Handling <dbl>,
#   Knotting_Binomial <chr>, Rotating_Duration <dbl>, Rotating_Prop <dbl>,
#   Rotating_Prop_Handling <dbl>, Ramming_Duration <dbl>, Ramming_Prop <dbl>,
#   Ramming_Prop_Handling <dbl>, Shaking_Duration <dbl>, Shaking_Prop <dbl>,
#   Shaking_Prop_Handling <dbl>, Crinkle_Duration <dbl>, Crinkle_Prop <dbl>,
#   Crinkle_Prop_Handling <dbl>, Tug_Duration <dbl>, Tug_Prop <dbl>, …
#Selecting the columns we want:
SelBehaviorDurationsP2 <- P2 %>%
  select(Treatment, Eel, Shaking_Duration, Ramming_Duration, Knotting_Duration, Rotating_Duration, Anchor_Duration, Crinkle_Duration) %>%
  drop_na()

head(SelBehaviorDurationsP2) 
# A tibble: 6 × 8
  Treatment Eel        Shaking_Duration Ramming_Duration Knotting_Duration
  <chr>     <chr>                 <dbl>            <dbl>             <dbl>
1 Small     Eelvis                0                    0                 0
2 Small     Snickers              1.13                 0                 0
3 Small     Bill                  0.504                0                 0
4 Small     Hamlet                0.403                0                 0
5 Small     Horatio               0                    0                 0
6 Small     Spaghettio            0                    0                 0
# ℹ 3 more variables: Rotating_Duration <dbl>, Anchor_Duration <dbl>,
#   Crinkle_Duration <dbl>
P2BehaviorDur_Long <- SelBehaviorDurationsP2 %>%
  pivot_longer(cols = c(Shaking_Duration, Ramming_Duration, Knotting_Duration, Rotating_Duration, Anchor_Duration, Crinkle_Duration), 
               names_to = "Behavior", 
               values_to = "Duration")

#Calculating Mean & SE - 0's included
P2BehaviorDurSE = P2BehaviorDur_Long %>%
  group_by(Treatment, Behavior) %>%
  summarize(mean = mean(Duration, na.rm = TRUE),
            sd = sd(Duration, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))
`summarise()` has grouped output by 'Treatment'. You can override using the
`.groups` argument.
P2BehaviorDurSE
# A tibble: 18 × 6
# Groups:   Treatment [3]
   Treatment Behavior             mean     sd     n     se
   <chr>     <chr>               <dbl>  <dbl> <int>  <dbl>
 1 Large     Anchor_Duration    2.22    4.14     20 0.925 
 2 Large     Crinkle_Duration   2.17    2.33     20 0.522 
 3 Large     Knotting_Duration  2.81    4.67     20 1.04  
 4 Large     Ramming_Duration   0       0        20 0     
 5 Large     Rotating_Duration  9.15    9.36     20 2.09  
 6 Large     Shaking_Duration   0.635   0.986    20 0.220 
 7 Medium    Anchor_Duration    2.34    2.95     20 0.660 
 8 Medium    Crinkle_Duration   6.31    7.00     20 1.57  
 9 Medium    Knotting_Duration  0       0        20 0     
10 Medium    Ramming_Duration   1.76    5.73     20 1.28  
11 Medium    Rotating_Duration 10.4    12.3      20 2.76  
12 Medium    Shaking_Duration   0.537   0.701    20 0.157 
13 Small     Anchor_Duration    6.15    6.23     17 1.51  
14 Small     Crinkle_Duration   9.40    9.14     17 2.22  
15 Small     Knotting_Duration  0       0        17 0     
16 Small     Ramming_Duration   0.0737  0.304    17 0.0737
17 Small     Rotating_Duration  7.47    8.41     17 2.04  
18 Small     Shaking_Duration   0.724   0.739    17 0.179 
#Making treatment a factor:
P2BehaviorDurSE$Treatment = factor(P2BehaviorDurSE$Treatment)

Stats - Duration of Behaviors

The next code is the statistical analysis, comparing the time spent engaged in each behavior (ramming, rotating, shaking, knotting, body anchoring [aka ‘crinkle’], and tail anchoring [aka ‘anchor’]) between the treatments (small, medium, large). To compare average duration spent engaged in each behavior, we fit the data for each behavior with a GLMM using a Tweedie distribution. The data is not normal and shows over disperson making other distributions less suitable. The Tweedie distribution encompasses a wide range of distribution shapes and is also well suited for modeling continuous positive data inclusive of 0 values.

Each behavior has an emoji header to quickly find each chunk.

#️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️
#                                                       BODY ANCHOR STATS 
#️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓⚓️⚓️⚓️⚓️️⚓⚓️️⚓️⚓️️⚓⚓️⚓️

# Pulling out the BODY ANCHORING duration data
BodyAnchorTweedie <- P2BehaviorDur_Long %>% 
  filter(Behavior == "Crinkle_Duration")


# Fit the model using glmmTMB
BodyAnchorTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                     data = BodyAnchorTweedie,
                                     family = tweedie)

# Print the model summary
summary(BodyAnchorTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: BodyAnchorTweedie

     AIC      BIC   logLik deviance df.resid 
   327.4    339.6   -157.7    315.4       51 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.144    0.3795  
Number of obs: 57, groups:  Eel, 8

Dispersion parameter for tweedie family (): 2.41 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)       0.7342     0.3159   2.324  0.02014 *  
TreatmentMedium   1.0391     0.3560   2.919  0.00351 ** 
TreatmentSmall    1.3985     0.3568   3.920 8.86e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
BodyAnchorDurEMMS <- emmeans(BodyAnchorTweedie_Model, "Treatment")
pairs(BodyAnchorDurEMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium   -1.039 0.356 Inf  -2.919  0.0098
 Large - Small    -1.399 0.357 Inf  -3.920  0.0003
 Medium - Small   -0.359 0.303 Inf  -1.188  0.4605

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(BodyAnchorDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Large      0.734 0.316 Inf     0.115      1.35  a    
 Medium     1.773 0.262 Inf     1.260      2.29   b   
 Small      2.133 0.265 Inf     1.613      2.65   b   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
#️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️
#                                                       TAIL ANCHOR STATS 
#️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓️⚓️⚓️⚓️⚓️️⚓⚓️⚓️⚓️⚓⚓⚓️⚓️⚓️⚓⚓️⚓️️⚓️⚓⚓️⚓️⚓️⚓️️


# Pulling out the TAIL ANCHOR duration data
TailAnchorTweedie <- P2BehaviorDur_Long %>% 
  filter(Behavior == "Anchor_Duration")

# Fit the model using glmmTMB
TailAnchorTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                    data = TailAnchorTweedie,
                                    family = tweedie)

# Print the model summary
summary(TailAnchorTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: TailAnchorTweedie

     AIC      BIC   logLik deviance df.resid 
   268.7    280.9   -128.3    256.7       51 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.4522   0.6725  
Number of obs: 57, groups:  Eel, 8

Dispersion parameter for tweedie family (): 3.08 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)  
(Intercept)       0.5565     0.4401   1.265   0.2060  
TreatmentMedium   0.1075     0.4505   0.239   0.8114  
TreatmentSmall    0.9096     0.4174   2.179   0.0293 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
TailAnchorDurEMMS <- emmeans(TailAnchorTweedie_Model, "Treatment")
pairs(TailAnchorDurEMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium   -0.107 0.450 Inf  -0.239  0.9691
 Large - Small    -0.910 0.417 Inf  -2.179  0.0748
 Medium - Small   -0.802 0.414 Inf  -1.936  0.1287

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(TailAnchorDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Large      0.557 0.440 Inf    -0.306      1.42  a    
 Medium     0.664 0.428 Inf    -0.174      1.50  a    
 Small      1.466 0.431 Inf     0.622      2.31  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# 😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫
#                                                             ROTATING STATS
# 😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫-😵‍💫


# Pulling out the ROTATING duration data
RotatingTweedie <- P2BehaviorDur_Long %>% 
  filter(Behavior == "Rotating_Duration")

# Fit the model using glmmTMB
RotatingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                    data = RotatingTweedie,
                                    family = tweedie)

# Print the model summary
summary(RotatingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: RotatingTweedie

     AIC      BIC   logLik deviance df.resid 
   386.5    398.8   -187.3    374.5       51 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.09566  0.3093  
Number of obs: 57, groups:  Eel, 8

Dispersion parameter for tweedie family (): 2.99 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)      2.17823    0.23995   9.078   <2e-16 ***
TreatmentMedium  0.08949    0.29390   0.304    0.761    
TreatmentSmall  -0.13163    0.32419  -0.406    0.685    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
RotatingDurEMMS <- emmeans(RotatingTweedie_Model, "Treatment")
pairs(RotatingDurEMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium  -0.0895 0.294 Inf  -0.304  0.9502
 Large - Small    0.1316 0.324 Inf   0.406  0.9132
 Medium - Small   0.2211 0.328 Inf   0.673  0.7791

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(RotatingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Small       2.05 0.269 Inf      1.52      2.57  a    
 Large       2.18 0.240 Inf      1.71      2.65  a    
 Medium      2.27 0.242 Inf      1.79      2.74  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# 🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏
#                                                                 RAMMING STATS
# 🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏🐏


# Pulling out the RAMMING duration data
RammingTweedie <- P2BehaviorDur_Long %>% 
  filter(Behavior == "Ramming_Duration")

# Fit the model using glmmTMB
RammingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                  data = RammingTweedie,
                                  family = tweedie)

# Print the model summary
summary(RammingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: RammingTweedie

     AIC      BIC   logLik deviance df.resid 
    59.5     71.8    -23.8     47.5       51 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 6.352    2.52    
Number of obs: 57, groups:  Eel, 8

Dispersion parameter for tweedie family (): 8.23 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)       -37.67   34505.99  -0.001    0.999
TreatmentMedium    36.22   34505.99   0.001    0.999
TreatmentSmall     32.94   34505.99   0.001    0.999
# Estimated marginal means and pairwise comparisons 
RammingDurEMMS <- emmeans(RammingTweedie_Model, "Treatment")
pairs(RammingDurEMMS)
 contrast       estimate       SE  df z.ratio p.value
 Large - Medium   -36.22 34505.99 Inf  -0.001  1.0000
 Large - Small    -32.94 34505.99 Inf  -0.001  1.0000
 Medium - Small     3.28     1.96 Inf   1.673  0.2158

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(RammingDurEMMS, Letters = letters)
 Treatment emmean       SE  df asymp.LCL asymp.UCL .group
 Large     -37.67 34505.99 Inf -67668.16  67592.83  a    
 Small      -4.73     3.68 Inf    -11.94      2.48  a    
 Medium     -1.44     3.27 Inf     -7.85      4.97  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
# 🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝
#                                                              SHAKING STATS 
# 🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝🤝

# Pulling out the SHAKING DURATION Data
ShakingTweedie <- P2BehaviorDur_Long %>% 
  filter(Behavior == "Shaking_Duration")


# Fit the model using glmmTMB
ShakingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                 data = ShakingTweedie,
                                 family = tweedie)

# Print the model summary
summary(ShakingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: ShakingTweedie

     AIC      BIC   logLik deviance df.resid 
   142.6    154.9    -65.3    130.6       51 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.2387   0.4886  
Number of obs: 57, groups:  Eel, 8

Dispersion parameter for tweedie family (): 0.798 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)  
(Intercept)      -0.5389     0.3122  -1.726   0.0843 .
TreatmentMedium  -0.2006     0.3642  -0.551   0.5817  
TreatmentSmall    0.1218     0.3559   0.342   0.7323  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Estimated marginal means and pairwise comparisons 
ShakingDurEMMS <- emmeans(ShakingTweedie_Model, "Treatment")
pairs(ShakingDurEMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium    0.201 0.364 Inf   0.551  0.8460
 Large - Small    -0.122 0.356 Inf  -0.342  0.9375
 Medium - Small   -0.322 0.365 Inf  -0.882  0.6515

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates 
emmeans:::cld.emmGrid(ShakingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Medium    -0.740 0.334 Inf     -1.39   -0.0858  a    
 Large     -0.539 0.312 Inf     -1.15    0.0729  a    
 Small     -0.417 0.322 Inf     -1.05    0.2131  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 3 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Graphics - Duration of Behaviors

# This graph shows the average DURATION spent in each behavior.
# 0's have been left in the data set 

#Naming the Future Facets:
Facet_Names_P2_BehaviorDur <- list(
  'Knotting_Duration' = "KNOTTING",
  'Ramming_Duration' = "RAMMING",
  'Rotating_Duration' = "ROTATING",
  'Shaking_Duration' = "SHAKING",
  'Crinkle_Duration' = "BODY ANCHOR",
  'Anchor_Duration' = "TAIL ANCHOR"
  )

Facet_labeller_BehaviorDur2<- function(variable,value){
  return(Facet_Names_P2_BehaviorDur[value])}


#Wrangling some pictures for the future
#Wrangling some pictures for the future
KnotPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Knotting.png") 
ShakePic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Shaking.png")
RotatePic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Rotating.png")
RamPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Ramming.png")
BodyAnchorPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/bodyanchor.png")
TailAnchorPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/tailanchor.png")


# Creating a data frame with pairwise significance 
P2_signif <- tribble(
  ~Treatment, ~Behavior,          ~Signif,
  "Large",    "Anchor_Duration", "b",    
  "Medium",   "Anchor_Duration", "b",   
  "Small",    "Anchor_Duration", "a",    
  "Large",    "Crinkle_Duration", "b",    
  "Medium",   "Crinkle_Duration", "a",   
  "Small",    "Crinkle_Duration", "a",    
  "Large",    "Knotting_Duration", "***",    
  "Medium",   "Knotting_Duration", "",   
  "Small",    "Knotting_Duration", "",    
  "Large",    "Ramming_Duration", "a",    
  "Medium",   "Ramming_Duration", "a",   
  "Small",    "Ramming_Duration", "a",   
  "Large",    "Rotating_Duration", "a",    
  "Medium",   "Rotating_Duration", "a",   
  "Small",    "Rotating_Duration", "a",    
  "Large",    "Shaking_Duration", "a",    
  "Medium",   "Shaking_Duration", "a",   
  "Small",    "Shaking_Duration", "a") %>% 
  left_join(transmute(P2BehaviorDurSE, Treatment, Behavior, mean, se),
            by = c("Treatment", "Behavior")) %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large"),
         Behavior = factor(Behavior),
         y = mean + se + 2)

# The graph:
P2_BehaviorDur = P2BehaviorDurSE %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           position = position_dodge()) +
  #Adding error bars: 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .5) + 
  #ADD SIGNIFICANCE TEXT
  geom_text(aes(y = y, label = Signif), P2_signif, size=6, face="bold") +
  #All other aesthetics:
  scale_fill_manual(values = c("darkslategray4", "darkslategray3", "azure2")) +
  scale_x_discrete(labels=c("Small", "Medium", "Large"))+
  ylim(0, 20)+
  labs(y = "Average Duration (Sec)", x = "Treatment") +
  facet_wrap(. ~ Behavior, labeller = Facet_labeller_BehaviorDur2) +
  theme_classic() +
  theme (axis.text.x = element_text(size = 15),
         axis.text.y =  element_text(size = 15),
         axis.title.x = element_text(size = 20, margin = margin(t = 10, r = 20, b = 0, l = 0)),
         axis.title.y = element_text(size = 20),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        plot.title = element_text(hjust =.5, size = 40, colour = "black"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face = "bold"))
Warning in geom_text(aes(y = y, label = Signif), P2_signif, size = 6, face =
"bold"): Ignoring unknown parameters: `face`
Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
arguments are now deprecated.
ℹ See labellers documentation.
P2_BehaviorDur

P2_BehaviorDur_Pics <- ggdraw(P2_BehaviorDur) +
draw_image(TailAnchorPic, x=0.31 , y=0.92, hjust = 1, vjust =1, height = .19, width
=.19) +
draw_image(BodyAnchorPic, x=0.62 , y=0.94, hjust = 1, vjust =1, height = .17, width
=.17) +
draw_image(KnotPic, x=0.95 , y=0.9, hjust = 1, vjust =1, height = .2, width
=.2 ) +
draw_image(RamPic, x=0.29 , y=0.45, hjust = 1, vjust =1, height = .17, width
=.17) +
draw_image(RotatePic, x=0.659 , y=0.54, hjust = 1, vjust =1, height = .23, width
=.23 ) +
draw_image(ShakePic, x=0.95 , y=0.45, hjust = 1, vjust =1, height = .2, width
=.2 )

P2_BehaviorDur_Pics 

#Saves to Reports folder within the project folder
ggsave(filename = "P2-Behavior Durations.png", plot = P2_BehaviorDur_Pics , width = 10, height = 7, dpi = 300)

—————————-

Experiment 2 - FINAL FIGURE

—————————

The code below combines the Trial Component Figure with the Behavior Duration Figure for Experiment 2

P2_FinalFigure <- ggarrange(P2_Dur_Graph, P2_BehaviorDur_Pics,
                    labels = c("A", "B"),
                    ncol = 1, nrow = 2,
                    heights = c(0.25, 0.5))
P2_FinalFigure

ggsave(filename = "P2_Final.png", plot = P2_FinalFigure, width = 11, height = 10, dpi = 300)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#~~~~~~~~~~~~~~~~~~~~ # BONUS FEATURES #~~~~~~~~~~~~~~~~~~~~

Experiment 1 & 2 Combined

The following code combines data from experiments 1 and 2. We wanted to compare the durations spent engaged in different behaviors between the open treatment (experiment 1) and all three diameter-based treatments from experiment 2 (small, medium, large)

Data Prep - Bonus

#Reading in the dataset

Bonus = read_csv(here::here("data/Phase2-SpatialContext+Open.csv")) 
Rows: 102 Columns: 30
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): Eel, Treatment, Knotting_Binomial
dbl (27): Treatment_Round, Trial_Duration, Handling_Duration, PT_Duration, P...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(Bonus) 
# A tibble: 6 × 30
  Eel     Treatment Treatment_Round Trial_Duration Handling_Duration PT_Duration
  <chr>   <chr>               <dbl>          <dbl>             <dbl>       <dbl>
1 Eelvis  Small                   1           50.4              49.3        1.17
2 Snicke… Small                   1           NA                NA         NA   
3 Bill    Small                   1           27.4              19.4        7.99
4 Hamlet  Small                   1           38.2              24.8       13.3 
5 Horatio Small                   1           30.4              17.9       12.5 
6 Spaghe… Small                   1           45.9              42.9        3.01
# ℹ 24 more variables: PT_Prop <dbl>, Knotting_Duration <dbl>,
#   Knotting_Prop <dbl>, Knotting_Prop_Handling <dbl>, Knotting_Binomial <chr>,
#   Rotating_Duration <dbl>, Rotating_Prop <dbl>, Rotating_Prop_Handling <dbl>,
#   Ramming_Duration <dbl>, Ramming_Prop <dbl>, Ramming_Prop_Handling <dbl>,
#   Shaking_Duration <dbl>, Shaking_Prop <dbl>, Shaking_Prop_Handling <dbl>,
#   Crinkle_Duration <dbl>, Crinkle_Prop <dbl>, Crinkle_Prop_Handling <dbl>,
#   Tug_Duration <dbl>, Tug_Prop <dbl>, Tug_Prop_Handling <dbl>, …
#Selecting the columns we want:
BonusSel <- Bonus %>%
  select(Treatment, Eel, Ramming_Duration, Shaking_Duration, Rotating_Duration, Knotting_Duration) %>%
  drop_na()

head(BonusSel) 
# A tibble: 6 × 6
  Treatment Eel        Ramming_Duration Shaking_Duration Rotating_Duration
  <chr>     <chr>                 <dbl>            <dbl>             <dbl>
1 Small     Eelvis                    0            0                14.1  
2 Small     Bill                      0            0.504             0.973
3 Small     Hamlet                    0            0.403             6.14 
4 Small     Horatio                   0            0                 0    
5 Small     Spaghettio                0            0                13.5  
6 Small     Jhonny                    0            0                 6.85 
# ℹ 1 more variable: Knotting_Duration <dbl>
#Making stuff a factor:
BonusSel$Treatment = factor(BonusSel$Treatment)
BonusSel$Eel = factor(BonusSel$Eel)


Bonus_Long <- BonusSel %>%
  pivot_longer(cols = c(Shaking_Duration, Ramming_Duration, Knotting_Duration, Rotating_Duration), 
               names_to = "Behavior", 
               values_to = "Duration")

head(Bonus_Long)
# A tibble: 6 × 4
  Treatment Eel    Behavior          Duration
  <fct>     <fct>  <chr>                <dbl>
1 Small     Eelvis Shaking_Duration     0    
2 Small     Eelvis Ramming_Duration     0    
3 Small     Eelvis Knotting_Duration    0    
4 Small     Eelvis Rotating_Duration   14.1  
5 Small     Bill   Shaking_Duration     0.504
6 Small     Bill   Ramming_Duration     0    
#Calculating Mean & SE 
BonusSE = Bonus_Long %>%
  group_by(Treatment, Behavior) %>%
  summarize(mean = mean(Duration, na.rm = TRUE),
            sd = sd(Duration, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))
`summarise()` has grouped output by 'Treatment'. You can override using the
`.groups` argument.
head(BonusSE)
# A tibble: 6 × 6
# Groups:   Treatment [2]
  Treatment Behavior           mean    sd     n    se
  <fct>     <chr>             <dbl> <dbl> <int> <dbl>
1 Large     Knotting_Duration 2.91  4.63     20 1.03 
2 Large     Ramming_Duration  0     0        20 0    
3 Large     Rotating_Duration 9.10  9.37     20 2.10 
4 Large     Shaking_Duration  0.666 0.974    20 0.218
5 Medium    Knotting_Duration 0     0        20 0    
6 Medium    Ramming_Duration  1.76  5.73     20 1.28

Bonus Stats

The following code compares the durations spent engaged in shaking, ramming, rotating, and knotting between the open and diameter based treatments (small, medium, large)

###################
## SHAKING DURATION
###################

BonusShaking <- Bonus_Long %>% 
  filter(Behavior == "Shaking_Duration") 


# Fit the model using glmmTMB
BonusShakingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                     data = BonusShaking,
                                     family = tweedie)

# Print the model summary
summary(BonusShakingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: BonusShaking

     AIC      BIC   logLik deviance df.resid 
   216.7    234.6   -101.3    202.7       89 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.3039   0.5512  
Number of obs: 96, groups:  Eel, 10

Dispersion parameter for tweedie family (): 0.844 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)  
(Intercept)      -0.5712     0.3241  -1.762    0.078 .
TreatmentMedium  -0.1975     0.3682  -0.536    0.592  
TreatmentOpen    -0.4890     0.3402  -1.438    0.151  
TreatmentSmall    0.2105     0.3768   0.559    0.576  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BonusShakingDurEMMS <- emmeans(BonusShakingTweedie_Model, "Treatment")
pairs(BonusShakingDurEMMS)
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium    0.198 0.368 Inf   0.536  0.9502
 Large - Open      0.489 0.340 Inf   1.438  0.4758
 Large - Small    -0.211 0.377 Inf  -0.559  0.9442
 Medium - Open     0.291 0.357 Inf   0.817  0.8467
 Medium - Small   -0.408 0.386 Inf  -1.057  0.7156
 Open - Small     -0.700 0.367 Inf  -1.905  0.2259

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
emmeans:::cld.emmGrid(BonusShakingDurEMMS, Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Open      -1.060 0.291 Inf     -1.63   -0.4888  a    
 Medium    -0.769 0.342 Inf     -1.44   -0.0991  a    
 Large     -0.571 0.324 Inf     -1.21    0.0640  a    
 Small     -0.361 0.339 Inf     -1.03    0.3041  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
###################
## RAMMING DURATION
###################

BonusRamming <- Bonus_Long  %>% 
  filter(Behavior == "Ramming_Duration")


BonusRammingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                     data = BonusRamming,
                                     family = tweedie)

# Print model summary
summary(BonusRammingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: BonusRamming

     AIC      BIC   logLik deviance df.resid 
   121.7    139.7    -53.9    107.7       89 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 7.929    2.816   
Number of obs: 96, groups:  Eel, 10

Dispersion parameter for tweedie family (): 5.93 

Conditional model:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)       -34.01   21206.20  -0.002    0.999
TreatmentMedium    31.56   21206.19   0.001    0.999
TreatmentOpen      32.09   21206.19   0.002    0.999
TreatmentSmall     28.39   21206.19   0.001    0.999
BonusRammingDurEMMS <- emmeans(BonusRammingTweedie_Model, "Treatment")
pairs(BonusRammingDurEMMS)
 contrast       estimate       SE  df z.ratio p.value
 Large - Medium  -31.558 21206.19 Inf  -0.001  1.0000
 Large - Open    -32.089 21206.19 Inf  -0.002  1.0000
 Large - Small   -28.387 21206.19 Inf  -0.001  1.0000
 Medium - Open    -0.531     0.96 Inf  -0.553  0.9457
 Medium - Small    3.171     1.81 Inf   1.754  0.2958
 Open - Small      3.702     1.98 Inf   1.874  0.2392

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
emmeans:::cld.emmGrid(BonusRammingDurEMMS, Letters = letters)
 Treatment emmean       SE  df asymp.LCL asymp.UCL .group
 Large     -34.01 21206.19 Inf -41597.38  41529.37  a    
 Small      -5.62     2.62 Inf    -10.75     -0.49  a    
 Medium     -2.45     1.95 Inf     -6.28      1.38  a    
 Open       -1.92     1.64 Inf     -5.12      1.29  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
####################
## ROTATING DURATION
####################

BonusRotating <-  Bonus_Long   %>% 
  filter(Behavior == "Rotating_Duration")

BonusRotatingTweedie_Model  <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                     data = BonusRotating,
                                     family = tweedie)
# Print model summary
summary(BonusRotatingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: BonusRotating

     AIC      BIC   logLik deviance df.resid 
   491.2    509.2   -238.6    477.2       89 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 0.1793   0.4235  
Number of obs: 96, groups:  Eel, 10

Dispersion parameter for tweedie family ():  2.3 

Conditional model:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      2.149971   0.242377   8.870  < 2e-16 ***
TreatmentMedium  0.056451   0.265611   0.213    0.832    
TreatmentOpen   -2.655312   0.335428  -7.916 2.45e-15 ***
TreatmentSmall   0.002388   0.294306   0.008    0.994    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Bonus_RotatingDurEMMS <- emmeans(BonusRotatingTweedie_Model, "Treatment")
pairs(Bonus_RotatingDurEMMS )
 contrast       estimate    SE  df z.ratio p.value
 Large - Medium -0.05645 0.266 Inf  -0.213  0.9966
 Large - Open    2.65531 0.335 Inf   7.916  <.0001
 Large - Small  -0.00239 0.294 Inf  -0.008  1.0000
 Medium - Open   2.71176 0.334 Inf   8.121  <.0001
 Medium - Small  0.05406 0.300 Inf   0.180  0.9979
 Open - Small   -2.65770 0.364 Inf  -7.293  <.0001

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
emmeans:::cld.emmGrid(Bonus_RotatingDurEMMS , Letters = letters)
 Treatment emmean    SE  df asymp.LCL asymp.UCL .group
 Open      -0.505 0.308 Inf     -1.11    0.0974  a    
 Large      2.150 0.242 Inf      1.67    2.6250   b   
 Small      2.152 0.272 Inf      1.62    2.6850   b   
 Medium     2.206 0.246 Inf      1.72    2.6895   b   

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
####################
## KNOTTING DURATION 
####################

Bonus_Knotting <- Bonus_Long %>% 
  filter(Behavior == "Knotting_Duration")

# Fit the gamma GLMM with repeated measures
BonusKnottingTweedie_Model <- glmmTMB(Duration ~ Treatment + (1 | Eel),
                                     data = Bonus_Knotting,
                                     family = tweedie)
# Print model summary
summary(BonusKnottingTweedie_Model)
 Family: tweedie  ( log )
Formula:          Duration ~ Treatment + (1 | Eel)
Data: Bonus_Knotting

     AIC      BIC   logLik deviance df.resid 
   227.3    245.3   -106.7    213.3       89 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 Eel    (Intercept) 1.795    1.34    
Number of obs: 96, groups:  Eel, 10

Dispersion parameter for tweedie family (): 5.84 

Conditional model:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)      9.464e-02  6.637e-01   0.143    0.887
TreatmentMedium -2.996e+01  1.759e+04  -0.002    0.999
TreatmentOpen    5.886e-01  4.870e-01   1.209    0.227
TreatmentSmall  -3.103e+01  3.173e+04  -0.001    0.999
Bonus_KnottingDurEMMS <- emmeans(BonusKnottingTweedie_Model, "Treatment")
pairs(Bonus_KnottingDurEMMS )
 contrast       estimate       SE  df z.ratio p.value
 Large - Medium   29.959 17587.11 Inf   0.002  1.0000
 Large - Open     -0.589     0.49 Inf  -1.209  0.6213
 Large - Small    31.030 31731.35 Inf   0.001  1.0000
 Medium - Open   -30.548 17587.11 Inf  -0.002  1.0000
 Medium - Small    1.071 36279.26 Inf   0.000  1.0000
 Open - Small     31.619 31731.35 Inf   0.001  1.0000

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
emmeans:::cld.emmGrid(Bonus_KnottingDurEMMS , Letters = letters)
 Treatment   emmean       SE  df asymp.LCL asymp.UCL .group
 Small     -30.9356 31731.35 Inf -62223.24   62161.4  a    
 Medium    -29.8644 17587.11 Inf -34499.96   34440.2  a    
 Large       0.0946     0.66 Inf     -1.21       1.4  a    
 Open        0.6833     0.57 Inf     -0.43       1.8  a    

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Bonus Graphics

The following code graphs the behaviors that showed interesting results (knotting and rotating)

BonusGraphData <- BonusSE %>% 
 filter(Behavior %in% c("Knotting_Duration", "Rotating_Duration"))

# Creating a data frame with pairwise significance 
Bonus_signif <- tribble(
  ~Treatment, ~Behavior,          ~Signif,
  "Open",     "Knotting_Duration", "a",
  "Large",    "Knotting_Duration", "a",    
  "Medium",   "Knotting_Duration",  "", 
  "Small",    "Knotting_Duration",  "",
  "Open",     "Rotating_Duration", "b",
  "Large",    "Rotating_Duration", "a",    
  "Medium",   "Rotating_Duration", "a",   
  "Small",    "Rotating_Duration", "a"
  ) %>% 
  left_join(transmute(BonusGraphData, Treatment, Behavior, mean, se),
            by = c("Treatment", "Behavior")) %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large", "Open"),
         Behavior = factor(Behavior),
         y = mean + se + 1)


#Naming the Future Facets:
Facet_Names_Bonus <- list(
  'Knotting_Duration' = "KNOTTING",
  'Rotating_Duration' = "ROTATING"
  )

Facet_labeller_Bonus <- function(variable,value){
  return(Facet_Names_Bonus[value])}


#Wrangling some pictures for the future
KnotPic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Knotting.png") 
RotatePic = here::here("/Users/mayamcelfish/Desktop/ArrRrrr/SpatialContext_PreyManipulation/images/Rotating.png")

# The graph:
BonusGraph = BonusGraphData %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large", "Open")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           position = position_dodge()) +
  #Adding error bars: 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .5) + 
  #ADD SIGNIFICANCE TEXT
  geom_text(aes(y = y, label = Signif), Bonus_signif, size=6, face="bold") +
  #All other aesthetics:
  scale_fill_manual(values = c("darkslategray4", "darkslategray3", "azure2", "white")) +
  scale_x_discrete(labels=c("Small", "Medium", "Large", "Open"))+
  ylim(0, 25)+
  labs(y = "Average Duration (Sec)", x = "Treatment") +
  facet_wrap(. ~ Behavior, labeller = Facet_labeller_Bonus) +
  theme_classic() +
  theme (axis.text.x = element_text(size = 15),
         axis.text.y =  element_text(size = 15),
         axis.title.x = element_text(size = 20, margin = margin(t = 10, r = 20, b = 0, l = 0)),
         axis.title.y = element_text(size = 20),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        plot.title = element_text(hjust =.5, size = 40, colour = "black"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face = "bold"))
Warning in geom_text(aes(y = y, label = Signif), Bonus_signif, size = 6, :
Ignoring unknown parameters: `face`
Warning: The `labeller` API has been updated. Labellers taking `variable` and `value`
arguments are now deprecated.
ℹ See labellers documentation.
BonusGraph

Bonus_Pics <- ggdraw(BonusGraph) +
draw_image(KnotPic, x=0.4 , y=0.86, hjust = 1, vjust =1, height = .25, width
=.25 ) +
draw_image(RotatePic, x=0.93 , y=0.9, hjust = 1, vjust =1, height = .3, width
=.3 ) 

Bonus_Pics

# This saves to the reports folder within the project
ggsave(filename = "BonusGraphic.png", plot = Bonus_Pics , width = 10, height = 5, dpi = 300)

🏎💨 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 🏎💨

ROTATION SPEED

🏎💨 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 🏎💨

Rotation Data Set-Up

The following code reads in the data set for rotating speed

#Reading in the dataset

Rot=read.csv(here::here("data/Rotation_Complete.csv"))

head(Rot) 
  Treatment Treatment_Exposure    Eel Rotation_Bout Rotation_Duration
1     Small                  1 Hamlet             1             3.250
2     Small                  1 Hamlet             2             1.739
3     Small                  1 Hamlet             3             2.239
4     Small                  2 Hamlet             1             1.688
5     Small                  2 Hamlet             2             1.003
6     Small                  2 Hamlet             3             1.502
  Opercular_Passes Opercular_Passes.2 Rotation_Second Notes
1                2                1.0       0.3076923    NA
2                3                1.5       0.8625647    NA
3                3                1.5       0.6699419    NA
4                2                1.0       0.5924171    NA
5                2                1.0       0.9970090    NA
6                3                1.5       0.9986684    NA

The followig code calculated the mean and standard error for rotating speed

#Making things factors:
Rot $ Treatment = factor(Rot $ Treatment)
Rot $ Eel = factor(Rot $ Eel)


#Calculating Mean & SE 
Rot_SE = Rot %>%
  group_by(Treatment) %>%
  summarize(mean = mean(Rotation_Second, na.rm = TRUE),
            sd = sd(Rotation_Second, na.rm = TRUE),
            n = n(),
            se = sd/sqrt(n))

Rot_SE
# A tibble: 4 × 5
  Treatment  mean    sd     n     se
  <fct>     <dbl> <dbl> <int>  <dbl>
1 Large     1.02  0.531    62 0.0674
2 Medium    0.679 0.365    58 0.0479
3 Open      1.42  0.434    51 0.0607
4 Small     0.692 0.350    58 0.0460

Rotation Speed Stats

To compare rotation speed between treatments, we log-transformed the data and employed a mixed-effects model with a gamma distribution. We incorporated individual eels as a repeated measure by including eel identity as a random effect to account for variability between individuals.

head(Rot)
  Treatment Treatment_Exposure    Eel Rotation_Bout Rotation_Duration
1     Small                  1 Hamlet             1             3.250
2     Small                  1 Hamlet             2             1.739
3     Small                  1 Hamlet             3             2.239
4     Small                  2 Hamlet             1             1.688
5     Small                  2 Hamlet             2             1.003
6     Small                  2 Hamlet             3             1.502
  Opercular_Passes Opercular_Passes.2 Rotation_Second Notes
1                2                1.0       0.3076923    NA
2                3                1.5       0.8625647    NA
3                3                1.5       0.6699419    NA
4                2                1.0       0.5924171    NA
5                2                1.0       0.9970090    NA
6                3                1.5       0.9986684    NA
# Log-transform the Rotation_Duration variable
Rot <- Rot %>% 
  mutate(Log_Rotation_Second= log(Rotation_Second + 1))  # Adding 1 to avoid log(0)


# Fit the gamma GLMM with repeated measures, adjusting optimization settings
RotGamma <- glmer(Rotation_Second ~ Treatment + (1 | Eel),
                  data = Rot,
                  family = Gamma(link = "log"))
 
# Print model summary
summary(RotGamma)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: Gamma  ( log )
Formula: Rotation_Second ~ Treatment + (1 | Eel)
   Data: Rot

     AIC      BIC   logLik deviance df.resid 
   188.7    208.4    -88.3    176.7      191 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8269 -0.6360 -0.1224  0.4818  4.3449 

Random effects:
 Groups   Name        Variance Std.Dev.
 Eel      (Intercept) 0.01339  0.1157  
 Residual             0.22742  0.4769  
Number of obs: 197, groups:  Eel, 10

Fixed effects:
                 Estimate Std. Error t value Pr(>|z|)    
(Intercept)     -0.004247   0.091480  -0.046    0.963    
TreatmentMedium -0.383990   0.089072  -4.311 1.63e-05 ***
TreatmentOpen    0.428373   0.110089   3.891 9.98e-05 ***
TreatmentSmall  -0.375937   0.095189  -3.949 7.84e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) TrtmnM TrtmnO
TreatmntMdm -0.496              
TreatmntOpn -0.440  0.414       
TretmntSmll -0.482  0.482  0.406
RotEMMS <- emmeans(RotGamma, "Treatment")
pairs(RotEMMS)
 contrast       estimate     SE  df z.ratio p.value
 Large - Medium  0.38399 0.0891 Inf   4.311  0.0001
 Large - Open   -0.42837 0.1101 Inf  -3.891  0.0006
 Large - Small   0.37594 0.0952 Inf   3.949  0.0005
 Medium - Open  -0.81236 0.1093 Inf  -7.434  <.0001
 Medium - Small -0.00805 0.0940 Inf  -0.086  0.9998
 Open - Small    0.80431 0.1125 Inf   7.148  <.0001

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 
emmeans:::cld.emmGrid(RotEMMS, Letters = letters)
 Treatment   emmean     SE  df asymp.LCL asymp.UCL .group
 Medium    -0.38824 0.0906 Inf    -0.566    -0.211  a    
 Small     -0.38018 0.0950 Inf    -0.566    -0.194  a    
 Large     -0.00425 0.0915 Inf    -0.184     0.175   b   
 Open       0.42413 0.1079 Inf     0.213     0.636    c  

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Graphics - Rotation Speeds

#Telling the bars to go from Small to Medium to Large instead of alphabetical 
Rot_SE$Treatment <- factor(Rot_SE$Treatment, levels=c("Small", "Medium", "Large","Open"))


# Step 1: create a data frame with pairwise significance 
# This is a dummy, adjust once you have the actual results
Rot_signif <- tribble(
  ~Treatment,   ~Signif,
  "Open",      "c",
  "Large",     "b",    
  "Medium",    "a",   
  "Small",     "a",    
 ) %>% 
  left_join(transmute(Rot_SE, Treatment, mean, se),
            by = c("Treatment")) %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large", "Open"),
         y = mean + se + .15)


Rot_Graph = Rot_SE %>% 
  mutate(Treatment = fct_relevel(Treatment, "Small", "Medium", "Large", "Open")) %>% 
  ggplot(aes(Treatment, mean, fill = Treatment)) + 
  geom_bar(stat ="identity", 
           color = "black", 
           width = 0.7,
           position = position_dodge()) +
  # error bar 
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), 
                width = .25) + 
  # STEP 2: ADD SIGNIFICANCE TEXT
  geom_text(aes(y = y, label = Signif), Rot_signif, size=6, face="bold") +
  scale_fill_manual(values = c("darkslategray4","darkslategray3", "azure2", "white")) +
  scale_x_discrete(labels=c("Small","Medium", "Large", "Open"))+
  labs(y = "Mean Rotation Speed (Rotations/Sec)", x = "Treatment", title = "") +
  # facet wrap + the label code from above for our new label names
  theme_classic() +
  theme (axis.text.x = element_text(size = 15),
         axis.text.y =  element_text(size = 15),
         axis.title.x = element_text(size = 20, margin = margin(t = 10, r = 20, b = 0, l = 0)),
         axis.title.y = element_text(size = 15),
        legend.position="none",
        panel.spacing.x = unit(2, "lines"),
        plot.title = element_text(hjust =.5, size = 40, colour = "black"),
        strip.background = element_rect(colour = "black", fill = "ivory3"),
        strip.placement = "inside",
        strip.text = element_text(size = 14, face = "bold"))
Warning in geom_text(aes(y = y, label = Signif), Rot_signif, size = 6, face =
"bold"): Ignoring unknown parameters: `face`
Rot_Graph 

# The graph saves to the Reports folder witn in the project folder
ggsave(filename = "Rotation Speed.png", plot = Rot_Graph, width = 5, height = 5, dpi = 300)