---
title: "Appendix S1. Computational details."
subtitle: "Modeling short-term energetic costs of sonar disturbance to cetaceans using high resolution foraging data"
author: 
 - "Max F. Czapanskiy"
 - "Matthew S. Savoca"
 - "William T. Gough"
 - "Paolo S. Segre"
 - "Danuta M. Wisniewska"
 - "David E. Cade"
 - "Jeremy A. Goldbogen"
date: "3/29/2021"
output: pdf_document
header-includes:
  - \usepackage[labelformat=empty]{caption}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(cowplot)
library(glue)
library(knitr)
library(kableExtra)
library(scales)
library(tidyverse)

set.seed(1651)
abbr_binom <- function(binom) {
  paste(str_sub(binom, 1, 1), 
        str_extract(binom, " .*"), 
        sep = ". ")
}
```

## Data

The sonar response energetic model uses three types of data: feeding rates, prey quality, and morphology. The feeding rate, morphology, and odontocete prey quality data were previously published by Goldbogen et al. (2019) [doi:10.1126/science.aax9044](https://doi.org/10.1126/science.aax9044). The rorqual prey quality data were published by Cade et al. (2021) [doi:10.1111/1365-2435.13763](https://doi.org/10.1111/1365-2435.13763). The reformatted data used in this analysis (prey_tbl.RDS, feeding_tbl.RDS, and morpho_tbl.RDS) are available through Dryad (see Data Availability statement).

```{r data}
prey <- readRDS("prey_tbl.RDS")
feeding <- readRDS("feeding_tbl.RDS")
morphology <- readRDS("morpho_tbl.RDS")

cbf_palette <- c(
  # blue
  "Delphinidae and Phocoenidae" = rgb(0, 114, 178, maxColorValue = 255),  
  # vermillion
  "Physeteridae and Ziphiidae" = rgb(213, 94, 0, maxColorValue = 255),
  # bluish green
  "Balaenopteridae" = rgb(0, 158, 115, maxColorValue = 255)
) 

morphology %>% 
  mutate(species = cell_spec(species, "latex", italic = TRUE)) %>% 
  select(Species = species, Family = family, Clade = clade, 
         `Length (m)` = length_m, `Mass (kg)` = mass_kg, 
         Abbreviation = abbr) %>% 
  kable("latex", booktabs = TRUE, escape = FALSE, 
        caption = "Table S1: Representative species morphology (reproduces Table 2 from main text).") %>%
  kable_styling(latex_options = c("scale_down"))
```

For the model, feeding rate ($r_f$) and energy acquired per prey capture event ($E_p$) distributions are empirical and lognormal, respectively. The following table combines the three data sources and provides quantile functions for $r_f$ and $E_p$ for each species.

```{r combine_data}
# Creates a quantile function based on the log-mean and log-sd of each species'
#   prey distribution
prey_fun <- prey %>% 
  mutate(q_Ep_kJ = map2(meanlnEp_lnkJ, 
                        sdlnEp_lnkJ, 
                        ~ function(p) qlnorm(p, .x, .y)))
# Creates an empirical quantile function for each species
feeding_fun <- feeding %>% 
  group_by(species) %>% 
  summarize(rf_h = list(rf_h)) %>% 
  mutate(q_rf_h = map(rf_h, ~ function(p) quantile(.x, p)))
# Combines the three data sources
species_data <- morphology %>% 
  left_join(prey_fun, by = "species") %>% 
  left_join(feeding_fun, by = "species")
```

## Undisturbed rate of energy acquisition

The undisturbed rate of energy acquisition ($P_a$) was modeled as the product of $r_f$ and $E_p$. The resulting distribution for blue whales (*Balaenoptera musculus*) is presented in the main manuscript (Fig. 2). The distributions for all species are below.

```{r Pa}
# Adapted from https://stackoverflow.com/a/45867076
label_sci <- function(breaks) {
  ifelse(
    breaks == 0, 
    "0", 
    parse(text = gsub("[+]", "", gsub("e", " %*% 10^", scientific(breaks))))
  )
}

# Plot the lognormal distribution of energy from prey
Ep_plot <- function(q_Ep_kJ, meanlnEp_lnkJ, sdlnEp_lnkJ) {
  tibble(x = q_Ep_kJ(c(0.001, 0.99))) %>% 
    ggplot(aes(x)) +
    stat_function(
      fun = dlnorm,
      args = list(meanlog = meanlnEp_lnkJ,
                  sdlog = sdlnEp_lnkJ),
      n = 1e3
    ) +
    expand_limits(x = 0) +
    scale_x_continuous(expression(italic(E[p]) ~ (kJ)), 
                       labels = label_sci) +
    theme_classic(base_size = 8) +
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.title.x = element_text(size = 8))
}

# Plot the empirical distribution of feeding rates
rf_plot <- function(q_rf_h) {
  tibble(x = q_rf_h(seq(0, 1, length.out = 1000))) %>% 
    ggplot(aes(x)) + 
    geom_histogram(bins = 20, 
                   boundary = 0, 
                   fill = "light gray", 
                   color = "black",
                   size = 0.2) +
    scale_x_continuous(expression(italic(r[f]) ~ ("hr"^{-1}))) +
    theme_classic(base_size = 8) +
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.title.x = element_text(size = 8))
}

# Plot the baseline rate of energy acquisition by sampling Ep and rf
Pa_plot <- function(species, q_Ep_kJ, q_rf_h, caption = TRUE) {
  Pa <- tibble(p_Ep = runif(1e3),
               p_rf = runif(1e3)) %>% 
    mutate(Ep = q_Ep_kJ(p_Ep),
           rf = q_rf_h(p_rf),
           Pa = Ep * rf)
  
  Pa_caption <- if (caption) {
    element_text(face = "italic")
  } else {
    element_blank()
  }
  
  Pa %>% 
    filter(Pa < quantile(Pa, 0.99)) %>% 
    ggplot(aes(Pa)) +
    geom_histogram(bins = 30,
                   boundary = 0, 
                   fill = "light gray", 
                   color = "black") +
    geom_vline(aes(xintercept = mean(Pa)),
               linetype = "dashed") +
    scale_x_continuous(expression(italic(P[a]) ~~ ("kJ" ~ "hr"^{-1})),
                       labels = label_sci) +
    labs(caption = species) +
    theme_classic(base_size = 12) +
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.title.x = element_text(size = 8),
          plot.caption = Pa_caption) 
}

# Combine plots (Ep and rf as insets in Pa)
arrange_Pa <- function(species, q_Ep_kJ, meanlnEp_lnkJ, sdlnEp_lnkJ, q_rf_h, 
                          caption = TRUE, ...) {
  Ep <- Ep_plot(q_Ep_kJ, meanlnEp_lnkJ, sdlnEp_lnkJ)
  rf <- rf_plot(q_rf_h)
  Pa <- Pa_plot(species, q_Ep_kJ, q_rf_h, caption)
  
  # Figure out where the insets should go
  scales <- layer_scales(Pa)
  xlim <- scales$x$get_limits()
  ylim <- scales$y$get_limits()
  grid_to_data <- function(xmin, xmax, ymin, ymax) {
    c(xmin = xlim[1] + xmin * (xlim[2] - xlim[1]),
      xmax = xlim[1] + xmax * (xlim[2] - xlim[1]),
      ymin = ylim[1] + ymin * (ylim[2] - ylim[1]),
      ymax = ylim[1] + ymax * (ylim[2] - ylim[1]))
  }
  Ep_coord <- grid_to_data(0.1, 0.535, 0.45, 0.95)
  rf_coord <- grid_to_data(0.555, 0.99, 0.45, 0.95)
  Pa +
    annotation_custom(ggplotGrob(Ep), 
                      xmin = Ep_coord["xmin"], xmax = Ep_coord["xmax"],
                      ymin = Ep_coord["ymin"], ymax = Ep_coord["ymax"]) +
    annotation_custom(ggplotGrob(rf), 
                      xmin = rf_coord["xmin"], xmax = rf_coord["xmax"],
                      ymin = rf_coord["ymin"], ymax = rf_coord["ymax"])
}

# Figure 2 of main text
fig2 <- species_data %>% 
  filter(abbr == "Bm") %>% 
  pmap(arrange_Pa, caption = FALSE) %>% 
  first()

Pa_plots <- species_data %>% 
  select(species, q_Ep_kJ, meanlnEp_lnkJ, sdlnEp_lnkJ, q_rf_h) %>% 
  pmap(arrange_Pa)
#walk(Pa_plots, print)
```

$r_f$ and mass-specific $E_p$ across species (Figure 3 in main text).

```{r fig3, fig.cap="Figure S1: Reproducing figure 3 from main text."}
fig3data <- species_data %>% 
  mutate(
    rf_mean = map_dbl(rf_h, mean),
    rf_1q = map_dbl(rf_h, quantile, p = 0.25),
    rf_3q = map_dbl(rf_h, quantile, p = 0.75),
    Ep_mean = exp(meanlnEp_lnkJ),
    Ep_1q = map2_dbl(meanlnEp_lnkJ, sdlnEp_lnkJ, ~ qlnorm(0.25, .x, .y)),
    Ep_3q = map2_dbl(meanlnEp_lnkJ, sdlnEp_lnkJ, ~ qlnorm(0.75, .x, .y)),
    Epmass_mean = Ep_mean / mass_kg,
    Epmass_1q = Ep_1q / mass_kg,
    Epmass_3q = Ep_3q / mass_kg,
    abbr = factor(abbr, levels = .$abbr)
  )
# Fig 3a
fig3a <- ggplot(fig3data, aes(abbr, rf_mean)) +
  geom_pointrange(aes(ymin = rf_1q, ymax = rf_3q, 
                      color = group, shape = group)) +
  scale_color_manual(values = cbf_palette) +
  labs(x = "", y = expression(italic(r[f])~~(hr^{-1}))) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")
# Fig 3b
fig3b <- ggplot(fig3data, aes(abbr, Epmass_mean)) +
  geom_pointrange(aes(ymin = Epmass_1q, ymax = Epmass_3q, 
                      color = group, shape = group)) +
  scale_color_manual(values = cbf_palette) +
  labs(x = "", y = expression("Mass-specific"~~italic(E[p])~~(kJ~kg^{-1}))) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")
# Fig 3
fig3 <- plot_grid(fig3a, fig3b, labels = "AUTO")
print(fig3)
```

```{r table5}
feeding_n <- feeding %>% 
  group_by(species) %>% 
  summarize(n_ind = n_distinct(id), 
            n_hrs = n(),
            .groups = "drop")
format_num <- function(x) {
  if (x < 1e3) {
    formatC(x, digits = 3, flag = "#")
  } else {
    exponent <- floor(log10(x))
    base <- signif(x / 10^exponent, 3)
    str_glue("{formatC(base, digits = 3, flag = \"#\")} $\\times 10^{exponent}$")
  }
}
format_Ep <- function(mean, q1, q3) {
  str_glue("{format_num(mean)} ({format_num(q1)} - {format_num(q3)})")
}
format_Pa <- function(q_Ep_kJ, q_rf_h) {
  p_Ep <- runif(1e3)
  p_rf <- runif(1e3) 
  Ep <- q_Ep_kJ(p_Ep)
  rf <- q_rf_h(p_rf)
  Pa <- Ep * rf
  Pa_mean <- format_num(mean(Pa))
  Pa_1q <- format_num(quantile(Pa, 0.25))
  Pa_3q <- format_num(quantile(Pa, 0.75))
  str_glue("{Pa_mean} ({Pa_1q} - {Pa_3q})")
}
format_Pa_m <- function(q_Ep_kJ, q_rf_h, mass_kg) {
  p_Ep <- runif(1e3)
  p_rf <- runif(1e3) 
  Ep <- q_Ep_kJ(p_Ep)
  rf <- q_rf_h(p_rf)
  Pa <- Ep * rf
  format_num(mean(Pa) / mass_kg)
}
table5 <- fig3data %>% 
  left_join(feeding_n, by = "species") %>% 
  transmute(species = str_glue("\\textit{{{abbr_binom(species)}}}\\\\({n_ind}, {n_hrs})"),
            rf = str_glue("{formatC(rf_mean, digits = 3, flag = \"#\")} ({rf_1q} - {rf_3q})"),
            Ep = pmap_chr(list(Ep_mean, Ep_1q, Ep_3q), format_Ep),
            Pa = map2_chr(q_Ep_kJ, q_rf_h, format_Pa),
            Pa_m = pmap_chr(list(q_Ep_kJ, q_rf_h, mass_kg), format_Pa_m))
colnames(table5) <- c(
  "Species",
  "$r_f$",
  "$E_p$",
  "$P_a$",
  "$P_a/m$"
)

table5 %>% 
  kable("latex", booktabs = TRUE, escape = FALSE,
        caption = "Table S2: Energy acquisition by species (reproduces Table 5 from main text).") %>% 
  kable_styling(latex_options = c("scale_down"))
```


## Energy expenditure due to elevated locomotor effort

Elevated locomotion is defined here as a speed increase ($U_f$) for a duration ($t_f$). The energetic cost associated with this behavior was modeled as the product of the increase in stroking frequency ($\Delta f$ in Hz), the mass-specific locomotor cost ($C_L$ in kJ stroke^-1^ kg^-1^), and the individual's mass. To incorporate uncertainty, both $\Delta f$ and $C_L$ were treated as gamma distributed variables with a mean equal to the best estimate and a shape parameter of 4. The following figure shows the estimated increase in mass-specific energy expenditure when $U_f$ is 5 m s^-1^. 

```{r deltaPe, fig.cap="Figure S2: Mass-specific change in energy expenditure."}
# Predict the change in fluking frequency (deltaff) and mass-specific locomotor
# costs (CL)
U_cruise_ms <- 1.5
U_flight_ms <- 5
# Change in fluking frequency in strokes per hour
ff_fun <- function(U, L, La = 0.2, St = 0.3) {
  St * U / L / La * 3600
}
# Mass-specific locomotor cost in kJ / kg / stroke
CL_fun <- function(m) (1.46 + 0.0005 * m) / 1000

# Estimate the quantile of mass-specific Pe (product of two gammas)
q_deltaPe <- function(p, deltaff_shape, deltaff_scale, CL_shape, CL_scale) {
  pmap_dbl(
    list(deltaff_shape, deltaff_scale, CL_shape, CL_scale),
    function(deltaff_shape, deltaff_scale, CL_shape, CL_scale) {
      tibble(
        deltaff = rgamma(1e6, shape = deltaff_shape, scale = deltaff_scale),
        CL = rgamma(1e6, shape = CL_shape, scale = CL_scale)
      ) %>% 
        mutate(deltaPe = deltaff * CL) %>% 
        pull(deltaPe) %>% 
        quantile(p)
    })
}

# Plot median and Q1-Q3 of mass-specific deltaPe
species_data %>% 
  mutate(
    ff_cruise = ff_fun(U_cruise_ms, length_m),
    ff_flight = ff_fun(U_flight_ms, length_m),
    mean_deltaff = ff_flight - ff_cruise,
    deltaff_shape = 4,
    deltaff_scale = mean_deltaff / deltaff_shape,
    mean_CL = CL_fun(mass_kg),
    CL_shape = 4,
    CL_scale = mean_CL / CL_shape,
    deltaPe_1q = q_deltaPe(0.25, deltaff_shape, deltaff_scale, CL_shape, CL_scale),
    deltaPe_med = q_deltaPe(0.5, deltaff_shape, deltaff_scale, CL_shape, CL_scale),
    deltaPe_3q = q_deltaPe(0.75, deltaff_shape, deltaff_scale, CL_shape, CL_scale),
    species = factor(species, levels = species[order(length_m)])
  ) %>% 
  ggplot(aes(x = species, y = deltaPe_med)) +
  geom_pointrange(aes(ymin = deltaPe_1q, ymax = deltaPe_3q)) +
  coord_flip() +
  labs(y = expression("Mass-specific " ~~ italic(Delta * P[e]) ~~ ("kJ" ~ "hr"^{-1} ~ "kg"^{-1}))) +
  theme_classic(base_size = 12) +
  theme(axis.text.y = element_text(face = "italic"),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_blank()) 
```

## Sensitivity analysis

We tested the model's sensitivity to energy acquisition ($E_p$, $r_f$) and expenditure ($\Delta f$, $C_L$) parameters in two behavioral scenarios. The scenarios were chosen to emphasize increased energy expenditure ($t_d$ = 1 hour, $t_f$ = 1 hour, $U_f$ = 5.0 m s^-1^) or lost consumption ($t_d$ = 4 hours, $t_f$ = 0.25 hours, $U_f$ = 3.0 m s^-1^). 

```{r sensitivity1}
# Energetic cost of sonar exposure
Esonar_fun <- function(td_hr, tf_hr, rf, Ep, delta_ff, CL, mass) {
  assimilation <- 0.84 # Lockyer 1981
  Pa <- assimilation * Ep * rf # kJ * hr-1 = kJ hr-1
  deltaPe <- delta_ff * CL * mass # strokes hr-1 * kJ kg-1 stroke-1 * kg = kJ hr-1
  Pa * td_hr + deltaPe * tf_hr # kJ hr-1 * hr + kJ hr-1 * hr = kJ
}
# Run a scenario
run_scenario <- function(species, scenario,
                         td_hr, tf_hr, Uf_ms,
                         q_rf_h, meanlnEp_lnkJ, sdlnEp_lnkJ, 
                         length_m, mass_kg, ...) {
  # When Uf equals cruising speed, delta_ff is 0, causing NaNs in qgamma
  qgamma2 <- function(p, shape, scale) {
    if (scale == 0)
      rep(0, length(p))
    else
      qgamma(p, shape = shape, scale = scale)
  }
  
  shape_arg <- 4
  delta_ff <- ff_fun(Uf_ms, length_m) - ff_fun(1.5, length_m)
  CL <- CL_fun(mass_kg)
  
  # Parameter list
  param <- c("rf_h", "Ep_kJ", "delta_ff", "CL")
  # Quantile functions
  q <- list(rf_h = q_rf_h,
            Ep_kJ = qlnorm,
            delta_ff = qgamma2,
            CL = qgamma)
  # Parameter arguments
  q_arg <- list(
    rf_h = list(),
    Ep_kJ = list(meanlog = meanlnEp_lnkJ, 
                 sdlog = sdlnEp_lnkJ),
    delta_ff = list(shape = shape_arg, scale = delta_ff / shape_arg),
    CL = list(shape = shape_arg, scale = CL / shape_arg)
  )
  # Generate Latin hypercube samples and run model
  pse::LHS(model = NULL, param, 5e2, q, q_arg, method = "random")$data %>% 
    mutate(
      esonar_kJ = Esonar_fun(td_hr, tf_hr, rf_h, Ep_kJ, delta_ff, CL, mass_kg)
    ) %>% 
    {cbind(tibble(species, scenario, td_hr, tf_hr, Uf_ms), .)}
}

# Run both behavioral scenarios
sensitivity_tbl <- 
  tribble(~scenario,     ~td_hr, ~tf_hr, ~Uf_ms,
          "flight",      1,      1.0,    5,
          "consumption", 4,      0.25,   3) %>% 
  crossing(species_data) %>% 
  pmap_dfr(run_scenario)
```

Model sensitivity to each parameter was quantified as the coefficient in the linear model $z(E) \sim z(r_f) + z(E_p) + z(\Delta f) + z(C_L)$ where $z()$ is the z-score. Results by guild are presented in the main text (Fig. 4). The following are the results for each species.

```{r sensitivity2, fig.cap="Figure S3: Model sensitivity by species (see Figure 4 in main text for sensitivity by guild.)"}
# Normalize parameters by z-score
zscore <- function(x) (x - mean(x)) / sd(x)
sensitivity_coef <- sensitivity_tbl %>% 
  group_by(species, scenario) %>% 
  mutate_at(vars(esonar_kJ, rf_h, Ep_kJ, delta_ff, CL), zscore) %>% 
  group_modify(function(data, keys) {
    # Fit linear model
    esonar_lm <- lm(esonar_kJ ~ rf_h + Ep_kJ + delta_ff + CL, data = data) 
    # Extract coefficients and confidence intervals
    esonar_coef <- coef(esonar_lm)
    esonar_ci <- as_tibble(confint(esonar_lm, level = 0.95),
                           rownames = "param")
    colnames(esonar_ci)[2:3] <- c("ci_min", "ci_max")
    
    cbind(esonar_coef, esonar_ci) %>% 
      select(param, esonar_coef, ci_min, ci_max) %>% 
      # drop intercept
      slice(-1)
  }) %>% 
  ungroup() %>% 
  mutate(
    param = factor(
      param, 
      levels = c("rf_h", "Ep_kJ", "delta_ff", "CL"),
      labels = c("r[f]", "E[p]", "Delta*f", "C[L]")
    ),
    species = factor(
      species,
      levels = species_data$species[order(species_data$length_m)]
    )
  )

# Plot sensitivity coefficients and confidence intervals for all species
param_lbls <- parse(text = levels(sensitivity_coef$param))
ggplot(sensitivity_coef, aes(species, esonar_coef, color = param)) +
  geom_pointrange(aes(ymin = ci_min, ymax = ci_max), 
                  position = position_dodge(1),
                  size = 0.25) +
  scale_color_viridis_d(labels = param_lbls) +
  coord_flip() +
  facet_wrap(~ scenario) +
  labs(y = "Sensitivity") +
  theme_classic(base_size = 12) +
  theme(axis.text.y = element_text(face = "italic"),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(hjust = 0),
        legend.title = element_blank(),
        strip.background = element_blank(),
        strip.text = element_blank()) 

# Sensitivity coefficients and CIs by guid
sensitivity_guilds <- sensitivity_tbl %>% 
  left_join(select(species_data, species, group), by = "species") %>% 
  mutate(group = factor(group, 
                        levels = c("Balaenopteridae",
                                   "Physeteridae and Ziphiidae",
                                   "Delphinidae and Phocoenidae"))) %>% 
  group_by(group, scenario) %>% 
  mutate_at(vars(esonar_kJ, rf_h, Ep_kJ, delta_ff, CL), zscore) %>% 
  group_modify(function(data, keys) {
    # Fit linear model
    esonar_lm <- lm(esonar_kJ ~ rf_h + Ep_kJ + delta_ff + CL, data = data) 
    # Extract coefficients and confidence intervals
    esonar_coef <- coef(esonar_lm)
    esonar_ci <- as_tibble(confint(esonar_lm, level = 0.95),
                           rownames = "param")
    colnames(esonar_ci)[2:3] <- c("ci_min", "ci_max")
    
    cbind(esonar_coef, esonar_ci) %>% 
      select(param, esonar_coef, ci_min, ci_max) %>% 
      # drop intercept
      slice(-1)
  }) %>% 
  ungroup() %>% 
  mutate(
    param = factor(
      param, 
      levels = c("rf_h", "Ep_kJ", "delta_ff", "CL"),
      labels = sprintf("italic(%s)", 
                       c("r[f]", "E[p]", "Delta*f", "C[L]"))
    )
  )
fig4 <- sensitivity_guilds %>% 
  group_by(scenario) %>% 
  group_map(
    function(data, ...) {
      ggplot(data,
             aes(x = param, y = esonar_coef, 
                 ymin = ci_min, ymax = ci_max, 
                 color = group)) +
        geom_pointrange(position = position_dodge(0.6), 
                        size = 1, 
                        fatten = 0.75) +
        coord_flip() +
        scale_x_discrete(
          labels = parse(text = levels(sensitivity_guilds$param))
        ) +
        scale_y_continuous(
          "Sensitivity",
          breaks = seq(0, 1, by = 0.2),
          labels = seq(0, 1, by = 0.2)
        ) +
        scale_color_manual(values = cbf_palette) +
        expand_limits(y = c(0, 1)) +
        theme_classic(base_size = 18) +
        theme(axis.title.y = element_blank(),
              legend.position = "none",
              strip.background = element_blank())
    }
  ) %>% 
  plot_grid(plotlist = ., labels = "AUTO")

# Table of sensitivity coefficients and CIs
coef_tbl <- sensitivity_coef %>% 
  mutate(abbr = abbr_binom(species),
         esonar_fmt = sprintf("%0.2f (%0.2f - %0.2f)", 
                              esonar_coef, ci_min, ci_max),
         species = as.character(species)) %>% 
  left_join(select(species_data, species, length_m), by = "species") %>% 
  arrange(length_m) %>% 
  select(scenario, abbr, param, esonar_fmt) %>% 
  pivot_wider(names_from = param, values_from = esonar_fmt) %>% 
  arrange(scenario) %>% 
  select(-scenario)
colnames(coef_tbl) <- c("", "$r_f$", "$E_p$", "$\\Delta f$", "$C_L$")

coef_tbl %>%
  kable("latex", booktabs = TRUE, escape = FALSE,
        caption = "Table S3: Model sensitivity by species.") %>% 
  kable_styling(latex_options = c("scale_down")) %>% 
  row_spec(0, align = "c") %>% 
  pack_rows("Lost consumption scenario", 1, 11) %>%
  pack_rows("Increased expenditure scenario", 12, 22)
```

## Relative energetic costs

Cross-species comparisons of energetic costs are complicated by the range of body sizes among cetaceans. A 10,000 kJ cost would have severe biological consequences for a harbor porpoise but not for a blue whale. Metabolic demands do not scale isometrically with body size, so the mass-specific energetic cost is not a proper comparison between species, either. We present a metric ($E^*$) that accounts for metabolic allometry by dividing the modeled energetic costs of four behavioral scenarios by the basal metabolic rate predicted by scaling relationships.

```{r estar}
# Calculate energy cost of four scenarios and E*
estar_scenarios <- tribble(
  ~scenario,   ~td_hr, ~tf_hr,  ~Uf_ms,
  "Mild",      1,      0.25,    2.5,
  "Strong",    2,      0.5,     3.5,
  "Extreme",   8,      2,       5
)

mr_tbl <- select(species_data, species, mass_kg, abbr, group) %>% 
  mutate(bmr_kleiber = 293.1 * mass_kg^0.75,
         # White and Seymour 2003, in ml O2 hr^-1 (use conversion factor of 
         # 20.1 kJ l O2^-1)
         bmr_white = 4.34 * (mass_kg * 1000)^0.67 / 1000 * 24 * 20.1,
         # Nagy et al 1999 in kJ day^-1
         fmr_nagy = 4.21 * (mass_kg * 1000)^0.772,
         # Maresh 2014 in kcal day^-1
         fmr_maresh = 732 * mass_kg^0.49 * 4.184) 

estar_tbl <- crossing(estar_scenarios, species_data) %>% 
  pmap_dfr(run_scenario) %>% 
  left_join(mr_tbl, 
            by = "species") %>% 
  mutate(estar = esonar_kJ / bmr_kleiber,
         estar_white = esonar_kJ / bmr_white,
         estar_nagy = esonar_kJ / fmr_nagy,
         estar_maresh = esonar_kJ / fmr_maresh,
         abbr = fct_reorder(abbr, mass_kg),
         scenario = factor(scenario, levels = c("Mild", "Strong", "Extreme")))

estar_summ <- estar_tbl %>% 
  group_by(group, species, abbr, mass_kg, scenario) %>% 
  summarize(estar_mean = mean(estar),
            estar_1q = quantile(estar, 0.25),
            estar_3q = quantile(estar, 0.75),
            estar_mean_white = mean(estar_white),
            estar_1q_white = quantile(estar_white, 0.25),
            estar_3q_white = quantile(estar_white, 0.75),
            estar_mean_nagy = mean(estar_nagy),
            estar_1q_nagy = quantile(estar_nagy, 0.25),
            estar_3q_nagy = quantile(estar_nagy, 0.75),
            estar_mean_maresh = mean(estar_maresh),
            estar_1q_maresh = quantile(estar_maresh, 0.25),
            estar_3q_maresh = quantile(estar_maresh, 0.75),
            .groups = "drop")
```

```{r estar_results, fig.cap="Figure S4: Relative immediate energetic costs of disturbance (reproduces Figure 5 from main text)."}
# Figure 5 from main text
fig5 <- ggplot(estar_summ, aes(mass_kg, estar_mean,
                       ymin = estar_1q, ymax = estar_3q,
                       color = group)) +
  geom_pointrange() +
  scale_color_manual(values = cbf_palette) +
  facet_grid(cols = vars(scenario)) +
  scale_x_continuous(trans = "log10",
                     breaks = trans_breaks("log10", function(x) 10^x),
                     labels = trans_format("log10", label_math(10^.x))) +
  scale_y_log10() +
  labs(x = expression(Mass ~ (kg)),
       y = expression(E * "*")) +
  theme_bw(base_size = 12) +
  theme(axis.text = element_text(size = 12),
        legend.position = "none", 
        strip.background = element_rect(color = "black", fill = "white"))
print(fig5)

# Behavioral scenarios for E* calculations
estar_scenarios %>% 
  set_names(c("Scenario", "$t_d~(hr)$", "$t_f~(hr)$", "$U_f~(m~s^{-1})$")) %>% 
  kable("latex", booktabs = TRUE, escape = FALSE,
        caption = "Table S4: Behavioral scenarios used to assess relative immediate energetic costs of disturbance.")

# Table of E* results
group_symbols <- tribble(
  ~group,                        ~group_sym,
  "Delphinidae and Phocoenidae", "\\ast",
  "Physeteridae and Ziphiidae",  "\\dag",
  "Balaenopteridae",             "\\ddag"
)
estar_table <- estar_summ %>%
  left_join(group_symbols, by = "group") %>% 
  mutate(
    estar_fmt = sprintf("%.3g (%.3g - %.3g)",
                        estar_mean, estar_1q, estar_3q),
    species_fmt = glue("$^{group_sym}\\textit{{{species}}}$")
  ) %>% 
  arrange(mass_kg) %>% 
  select(species_fmt, scenario, estar_fmt) %>% 
  pivot_wider(names_from = scenario, values_from = estar_fmt)
colnames(estar_table) <- c(
  "Species", "Mild", "Strong", "Extreme"
)
kable(estar_table, "latex", booktabs = TRUE, escape = FALSE,
        caption = "Table S5: Relative immediate energetic costs of disturbance by species and scenario.") %>%
  kable_styling(latex_options = c("scale_down")) %>% 
  add_header_above(c(" " = 1, "$E^*$ by behavioral scenario" = 3),
                   escape = FALSE) %>% 
  add_footnote(group_symbols$group, notation = "symbol") %>% 
  collapse_rows(columns = 1, latex_hline = "major", valign = "top")
```

```{r applications}
tribble(
  ~species,                 ~scenario, ~td_hr, ~tf_hr, ~Uf_ms,
  "Megaptera novaeangliae", "HI",      0.0,    0.3,    2.8,
  "Balaenoptera musculus",  "CA",      0.25,   0.2,    3.0,
  "Ziphius cavirostris",    "CA",      7.6,    1.7,    3.1
) %>% 
  left_join(species_data, by = "species") %>% 
  pmap_dfr(run_scenario) %>% 
  group_by(species, scenario) %>% 
  summarize(E_kJ_med = median(esonar_kJ),
            E_kJ_1Q = quantile(esonar_kJ, 0.25),
            E_kJ_3Q = quantile(esonar_kJ, 0.75),
            .groups = "drop") %>% 
  mutate_if(is.numeric, ~ format(., digits = 3, scientific = TRUE)) %>% 
  mutate(E = glue("{E_kJ_med} ({E_kJ_1Q} - {E_kJ_3Q})")) %>% 
  select(Species = species, Location = scenario, E) %>% 
  slice(c(2, 1, 3))
```

```{r alternative_estar, fig.cap=c("Figure S5: Relative immediate energetic costs of disturbance relative to basal metabolic rate (as estimated by allometric equation in (White and Seymour, 2003)).", "Figure S6: Relative immediate energetic costs of disturbance relative to field metabolic rate (as estimated by allometric equation for eutherians in (Nagy et al., 1999)).", "Figure S7: Relative immediate energetic costs of disturbance relative to marine mammal field metabolic rate (as estimated by allometric equation in (Maresh 2014)).")}
# Reproduce figure 5 from main text using alternative metabolic rate metrics
fig5_white <- ggplot(estar_summ, aes(mass_kg, estar_mean_white,
                                     ymin = estar_1q_white, 
                                     ymax = estar_3q_white,
                                     color = group)) +
  geom_pointrange() +
  scale_color_manual(values = cbf_palette) +
  facet_grid(cols = vars(scenario)) +
  scale_x_continuous(trans = "log10",
                     breaks = trans_breaks("log10", function(x) 10^x),
                     labels = trans_format("log10", label_math(10^.x))) +
  scale_y_log10() +
  labs(x = expression(Mass ~ (kg)),
       y = expression(E * "*")) +
  theme_bw(base_size = 12) +
  theme(axis.text = element_text(size = 12),
        legend.position = "none",
        strip.background = element_rect(color = "black", fill = "white"))
print(fig5_white)

fig5_nagy <- ggplot(estar_summ, aes(mass_kg, estar_mean_nagy,
                                     ymin = estar_1q_nagy, 
                                     ymax = estar_3q_nagy,
                                     color = group)) +
  geom_pointrange() +
  scale_color_manual(values = cbf_palette) +
  facet_grid(cols = vars(scenario)) +
  scale_x_continuous(trans = "log10",
                     breaks = trans_breaks("log10", function(x) 10^x),
                     labels = trans_format("log10", label_math(10^.x))) +
  scale_y_log10() +
  labs(x = expression(Mass ~ (kg)),
       y = expression(E * "*")) +
  theme_bw(base_size = 12) +
  theme(axis.text = element_text(size = 12),
        legend.position = "none",
        strip.background = element_rect(color = "black", fill = "white"))
print(fig5_nagy)

fig5_maresh <- ggplot(estar_summ, aes(mass_kg, estar_mean_maresh,
                                     ymin = estar_1q_maresh, 
                                     ymax = estar_3q_maresh,
                                     color = group)) +
  geom_pointrange() +
  scale_color_manual(values = cbf_palette) +
  facet_grid(cols = vars(scenario)) +
  scale_x_continuous(trans = "log10",
                     breaks = trans_breaks("log10", function(x) 10^x),
                     labels = trans_format("log10", label_math(10^.x))) +
  scale_y_log10() +
  labs(x = expression(Mass ~ (kg)),
       y = expression(E * "*")) +
  theme_bw(base_size = 12) +
  theme(axis.text = element_text(size = 12),
        legend.position = "none",
        strip.background = element_rect(color = "black", fill = "white"))
print(fig5_maresh)
```

```{r session_info}
sessionInfo()
```


