Skip to content

Instantly share code, notes, and snippets.

@FlukeAndFeather

FlukeAndFeather/sonar_supplement.Rmd

Last active February 22, 2021 23:24
Show Gist options
  • Save FlukeAndFeather/6c0ae9a42558d5e899354a8701db8f4a to your computer and use it in GitHub Desktop.
Save FlukeAndFeather/6c0ae9a42558d5e899354a8701db8f4a to your computer and use it in GitHub Desktop.
Download ZIP
Supporting Information for "Large baleen and small toothed whales face greatest energetic consequences from sonar disturbance"
Raw
sonar_supplement.Rmd
---
title: "Supporting Information"
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: "2/22/2021"
output: pdf_document
---
```{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).
```{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) %>%
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}
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) %>%
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}
# 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}
# 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) %>%
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}
# 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)
# 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) %>%
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}
# 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 * "*"),
caption = "E* relative to BMR (White and Seymour 2003)") +
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 * "*"),
caption = "E* relative to FMR (Nagy et al. 1999)") +
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 * "*"),
caption = "E* relative to FMR (Maresh 2014)") +
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)
```
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment