The objective with this analysis was to explore the allometric relationship between breathing frequency and body mass for mammals that live in different habitats.
Unlike other studies, We aimed to only use data that were sampled under basal conditions; 1) adult animals, 2) not pregnant or lactating, 3) post-prandial (fasted), 4) at rest, and 5) in their thermoneutral zone. Although, we included animals at different activity levels, the final analysis focused on those that were inactive, so activity level 1. Also, not all animals may have fasted long enough to be post-prandial, like ruminants.
gf_point(log_fr ~ log_Mb | habitat, alpha=0.1, data = breath_data)
gf_point(log_fr ~ log_Mb, data = breath_data, color=~habitat)|>gf_labs(x ="Log 10 Body mass (kg)",y ="Log 10 Breaths per Minute",color ="Habitat")
Model Fitting
This model that includes consideration (phylogenetic contrast) of order, genus and species via a set of nested random intercepts, plus a separate random intercept for the individual animal. The model includes body mass, activity, temperature of either water or air (where the animal was measured) and habitat. The model also accounts for differences in slopes (for the body mass effect) for animals of different habitats.
These graphical checks were carried out to confirm that model conditions (linearity, plus constant variance, independence, and normality of residuals) were met.
The ANOVA output shows a strong evidence of association between breath rate and body mass, habitat, and activity level, and provides moderate evidence of an interaction between habitat and body mass. The data also provides strong evidence of a difference by location, but note that only semi-aquatic mammals are measured in both locations (as all aquatic are only measured in water and all terrestrial in air). So, part of the location effect may be driven by overall terrestrial-aquatic differences. Surprisingly, location suggests that breathing frequency is lower in water than on land for semi-aquatic mammals.
To dig deeper into the relationship between breath rate, mass, and habitat, we used the R package emmeans and its function emtrends() to do a post-hoc comparison of the slopes for each habitat.
For each pair of habitats, we will test the null hypothesis that the difference in slopes is 0.
(So, a small p-value for a pair means the slopes of those two habitats are different.)
Since we are doing several hypothesis tests together, the chances of a false positive are elevated; to compensate, emtrends() has used Tukey’s method to adjust the p-values before reporting.
plot(slope_comparison$emtrends) +labs(x ="log(Mass (kg)) Trend", y ="Habitat")
The top part of the output gives a slope for each habitat with CIs. The lower “contrasts” part is where we’ll focus now: it gives the p-values for the tests comparing each pair of habitats.
We have moderate evidence that terrestrial habitats are different from the others. But we have no evidence that the slopes are different between aquatic and semiaquatic.
Model Predictions
Predictions by Habitat (and Location)
For the prediction plots shown below, the other predictors in the model but not shown in the plot were held fixed at the following values:
activity level: 1
location: water for aquatic, air for terrestrial, both for semiaquatic
We also add in data from previous studies by Mortola and He for comparison.
Open squares are data points from He et al. Open triangles are data points from Mortola et al. Grey data points are taken in water, and black ones in air. Lines and ribbons show model predictions, with the solid lines indicating measurements on land and dotted lines in water.
Predictions by Habitat (and Location) - Interactive Version
breath_data <- breath_data |>mutate(individual_info =paste0(individual_name, ' (', common_name, ', ', genus, ' ', species, ')'))gf_point(breath_rate ~ body_mass | habitat, data = breath_data,alpha =0.3, label =~individual_info, color =~order) |># with log scales on both axesgf_refine(scale_x_log10(), scale_y_log10()) |>gf_labs(x ='Body Mass (kg)',y ='Breaths per Minute') |># add the predictions nowgf_line(breaths ~ mass | habitat,linetype =~location,data = natural_preds,inherit =FALSE) |>gf_ribbon(conf_low + conf_hi ~ mass | habitat,linetype =~location,data = natural_preds,inherit =FALSE) |># add in He species -- grey squaresgf_point(breath_rate ~ body_mass | habitat,inherit =FALSE,data = he,color ='grey44',label =~common_name,shape =15) |># add in Mortola species -- grey trianglesgf_point(breath_rate ~ body_mass | habitat,inherit =FALSE,data = mortola,color ='grey44',label =~common_name,shape =17) |>gf_refine(guides(linetype ="none")) |>ggplotly()
Source Code
---title: "Allometry of Mammalian Breathing Rates"author: "Elliot Stielstra, Ethan Wilstermann, Simon Rylaarsdam, Stacy DeRuiter, Andreas Fahlman"date: "last-modified"format: html: embed-resources: true code-tools: true toc: trueeditor: source---```{r}#| label: setup #| include: falselibrary(tidyverse)library(ggformula)library(glmmTMB)library(ggeffects)library(readxl)library(plotly)library(DHARMa)library(emmeans)knitr::opts_chunk$set(# folder where images will be savedfig.path ='figures/breath-rate-',# file type for image files (can also be eps, jpeg, ...)dev =c('png', 'tiff', 'postscript'),# resolution for image filesdpi =300,warning =FALSE,message =FALSE,fig.width =6.5, fig.height =3.5)theme_set(theme_minimal(base_size =12))```# Data AccessRaw data files can be found at:<https://stacyderuiter.github.io/breath-allometry/data/breath-rate-FINAL-clean.csv><https://stacyderuiter.github.io/breath-allometry/data/additional-species.xlsx>The public github repository containing this work is at:<https://github.com/stacyderuiter/breath-allometry># Analysis## Context and GoalsThe objective with this analysis was to explore the allometric relationship between breathing frequency and body mass for mammals that live in different habitats. Unlike other studies, We aimed to only use data that were sampled under basal conditions; 1) adult animals, 2) not pregnant or lactating, 3) post-prandial (fasted), 4) at rest, and 5) in their thermoneutral zone. Although, we included animals at different activity levels, the final analysis focused on those that were inactive, so activity level 1. Also, not all animals may have fasted long enough to be post-prandial, like ruminants. ## Data Preparation```{r}breath_data <-read_csv('data/breath-rate-FINAL.csv',show_col_types =FALSE,trim_ws =TRUE,locale = readr::locale(encoding ="latin1")) #extract the variables of interestbreath_data <- breath_data |>select(log_fr, log_Mb, habitat, order, family, genus, species, breath_rate, body_mass, sex, activity_level, individual_name, common_name, age, temp, location,`Thermoneutral range`, `Normal temperature range`) |>mutate(across(c(order, family, genus, individual_name), function(x) str_replace_all(x, "[^[:alnum:]]", ""))) |>mutate(across(c(species, common_name), function(x) str_replace_all(x, "[^[:alnum:]]", " "))) |>mutate(common_name =ifelse(species =="orca", "Killer Whale", common_name)) |>mutate(across(c(individual_name, common_name, habitat, location), str_to_title)) |>mutate(common_name =str_replace_all(common_name, pattern ="Harbour", replacement ="Harbor")) |>mutate(common_name =str_replace_all(common_name, pattern ="Sabel", replacement ="Sable")) |>mutate(common_name =str_replace_all(common_name, pattern ="Giraff", replacement ="Giraffe")) |>mutate(common_name =str_replace_all(common_name, pattern ="Giraffee", replacement ="Giraffe")) |>mutate(common_name =str_replace_all(common_name, pattern ="Asiatisk", replacement ="Asian")) |>mutate(common_name =ifelse(common_name =="Brown Fur Seal", "South African Fur Seal", common_name)) |>mutate(activity_level =factor(activity_level),habitat =factor(habitat),location =factor(location)) |>drop_na(log_fr, log_Mb, habitat, activity_level, temp, location, order, family, genus, species, individual_name)breath_data <- breath_data |>arrange(order, common_name)write_csv(breath_data, 'data/breath-rate-FINAL-clean.csv')glimpse(breath_data)```## Data Exploration```{r, simple-data-scatterplot}gf_point(log_fr ~ log_Mb | habitat, alpha=0.1, data = breath_data) ``````{r, scatter-color-by-habitat}gf_point(log_fr ~ log_Mb, data = breath_data, color= ~habitat)|> gf_labs(x = "Log 10 Body mass (kg)", y = "Log 10 Breaths per Minute", color = "Habitat")```## Model FittingThis model that includes consideration (phylogenetic contrast) of order, genus and species via a set of nested random intercepts, plus a separate random intercept for the individual animal. The model includes body mass, activity, temperature of either water or air (where the animal was measured) and habitat. The model also accounts for differences in slopes (for the body mass effect) for animals of different habitats.```{r}breath_model <-glmmTMB(log_fr ~ log_Mb * habitat + activity_level + temp + location + (1| order/family/genus/species) + (1| individual_name), data = breath_data)summary(breath_model)```## Model AssessmentThese graphical checks were carried out to confirm that model conditions (linearity, plus constant variance, independence, and normality of residuals) were met.```{r}breath_data <- breath_data |>mutate(preds =predict(breath_model),resids =resid(breath_model))``````{r, resid-vs-fitted}gf_point(resids ~ preds, data = breath_data, alpha = 0.2)``````{r, scaled-resid}sim_res <- simulateResiduals(breath_model)plotResiduals(sim_res, quantreg = FALSE)``````{r, resid-acf}acf(resid(breath_model))``````{r, resid-hist}gf_histogram(~resids, data = breath_data, bins = 20)```## Inference```{r}car::Anova(breath_model)```The ANOVA output shows a strong evidence of association between breath rate and body mass, habitat, and activity level, and provides moderate evidence of an interaction between habitat and body mass. The data also provides strong evidence of a difference by location, but note that only semi-aquatic mammals are measured in both locations (as all aquatic are only measured in water and all terrestrial in air). So, part of the location effect may be driven by overall terrestrial-aquatic differences. Surprisingly, location suggests that breathing frequency is lower in water than on land for semi-aquatic mammals. To dig deeper into the relationship between breath rate, mass, and habitat, we used the R package `emmeans` and its function `emtrends()` to do a post-hoc comparison of the slopes for each habitat. For each pair of habitats, we will test the null hypothesis that the *difference in slopes* is 0. (So, a small p-value for a pair means the slopes of those two habitats are different.) Since we are doing several hypothesis tests together, the chances of a false positive are elevated; to compensate, `emtrends()` has used Tukey's method to adjust the p-values before reporting.```{r}slope_comparison <-emtrends(breath_model, pairwise ~ habitat,var ="log_Mb")slope_comparison``````{r}slope_comparison$emtrends |>data.frame() |>rename(Habitat = habitat,`log(Mass (kg)) Trend`=`log_Mb.trend`) |>mutate(across(where(is.numeric), function(x) round(x, digits =3))) |>mutate(`95% CI`=paste0("(", lower.CL, ", ", upper.CL, ")")) |>select(Habitat, `log(Mass (kg)) Trend`, SE, df, `95% CI`) |> flextable::flextable()plot(slope_comparison$emtrends) +labs(x ="log(Mass (kg)) Trend", y ="Habitat")```The top part of the output gives a slope for each habitat with CIs. The lower "contrasts" part is where we'll focus now: it gives the p-values for the tests comparing each pair of habitats.We have moderate evidence that terrestrial habitats are different from the others. But we have no evidence that the slopes are different between aquatic and semiaquatic.## Model Predictions### Predictions by Habitat (and Location)For the prediction plots shown below, the other predictors in the model but not shown in the plot were held fixed at the following values:- activity level: 1- location: water for aquatic, air for terrestrial, both for semiaquaticWe also add in data from previous studies by Mortola and He for comparison.```{r}mortola <-read_excel('data/additional-species.xlsx',sheet ='mortola') |>mutate(log_Mb =log10(`Body mass`),body_mass =`Body mass`,log_fr =log10(fr),breath_rate = fr,common_name = species,across(c(common_name, habitat), str_to_title)) he <-read_excel('data/additional-species.xlsx',sheet ='He') |>rename(log_Mb =`Log of Mass`,body_mass =`Mass (kg)`,breath_rate =`Breathing Frequency (br/min)` ) |>mutate(log_fr =log10(breath_rate),common_name =`Common Name`,across(c(common_name, habitat), str_to_title))``````{r, data-plus-model-static-with-n}colrs <- c('black', 'grey70')log_preds <- ggpredict(breath_model, terms = c('log_Mb', 'habitat', 'location'), condition = c(activity_level = 1))log_preds <- log_preds |> filter((facet == 'Air' & group %in% c('Terrestrial', 'Semiaquatic')) | (facet == 'Water' & group %in% c('Aquatic', 'Semiaquatic')))individs_by_habitat <- breath_data |> select(log_fr, log_Mb, habitat, activity_level, temp, location, order, family, genus, species, individual_name,) |> drop_na() |> group_by(order, family, genus, species, individual_name, habitat) |> summarize(individual_name = first(individual_name)) |> ungroup() |> group_by(habitat) |> summarize(n_individs = n())measurements_by_habitat <- breath_data |> select(log_fr, log_Mb, habitat, activity_level, temp, location, order, family, genus, species, individual_name,) |> drop_na() |> group_by(habitat) |> summarize(n_measurements = n())habitat_labs <- paste0(levels(breath_data$habitat), '\n(', measurements_by_habitat$n_measurements, ' points, ', individs_by_habitat$n_individs, ' animals)')names(habitat_labs) <- levels(breath_data$habitat)natural_preds <- data.frame(mass = 10^log_preds$x, breaths = 10^log_preds$predicted, conf_low = 10^log_preds$conf.low, conf_hi = 10^log_preds$conf.high, habitat = log_preds$group, location = log_preds$facet )# plot datagf_point(breath_rate ~ body_mass, color = ~location, data = breath_data, alpha = 0.6 ) |> gf_facet_grid(~habitat , labeller = labeller(habitat = habitat_labs) ) |> # with log scales on both axes gf_refine(scale_x_log10(), scale_y_log10(), scale_color_manual("", values = colrs)) |> gf_labs(x = 'Body Mass (kg)', y = 'Breaths per Minute') |> # add in He species -- open squares gf_point(breath_rate ~ body_mass | habitat, # inherit = FALSE, data = he, color = 'black', label = ~common_name, shape = 0) |> # 15) |> # add in Mortola species -- open triangles gf_point(breath_rate ~ body_mass | habitat, # inherit = FALSE, data = mortola, color = 'black', label = ~common_name, shape = 2) |> #17) |> # add the predictions now gf_line(breaths ~ mass, linetype = ~location, data = natural_preds, inherit = FALSE) |> gf_facet_grid(~habitat , labeller = labeller(habitat = habitat_labs) ) |> gf_ribbon(conf_low + conf_hi ~ mass | habitat, linetype = ~location, # fill = ~location, data = natural_preds, inherit = FALSE, alpha = 0.2) |> gf_facet_grid(~habitat , labeller = labeller(habitat = habitat_labs) ) |> gf_refine(guides(color = 'none', linetype = 'none'))```Open squares are data points from He et al. Open triangles are data points from Mortola et al. Grey data points are taken in water, and black ones in air. Lines and ribbons show model predictions, with the solid lines indicating measurements on land and dotted lines in water.### Predictions by Habitat (and Location) - Interactive Version```{r, data-plus-model-interactive}breath_data <- breath_data |> mutate(individual_info = paste0(individual_name, ' (', common_name, ', ', genus, ' ', species, ')'))gf_point(breath_rate ~ body_mass | habitat, data = breath_data, alpha = 0.3, label = ~individual_info, color = ~order) |> # with log scales on both axes gf_refine(scale_x_log10(), scale_y_log10()) |> gf_labs(x = 'Body Mass (kg)', y = 'Breaths per Minute') |> # add the predictions now gf_line(breaths ~ mass | habitat, linetype = ~location, data = natural_preds, inherit = FALSE) |> gf_ribbon(conf_low + conf_hi ~ mass | habitat, linetype = ~location, data = natural_preds, inherit = FALSE) |> # add in He species -- grey squares gf_point(breath_rate ~ body_mass | habitat, inherit = FALSE, data = he, color = 'grey44', label = ~common_name, shape = 15) |> # add in Mortola species -- grey triangles gf_point(breath_rate ~ body_mass | habitat, inherit = FALSE, data = mortola, color = 'grey44', label = ~common_name, shape = 17) |> gf_refine(guides(linetype = "none")) |> ggplotly()```
