In robschick/tangled: Creates data objects and plots for Knowlton, Schick et al. manuscript on entanglement

The goal of this vignette is to describe median health during the (shortened) entanglement windows as a function of gear-carrying status and injury severity. To do this, we will create a series of paired boxplots showing health anomaly during the entanglment window for both non-reproductive animals and reproductive females. This will also include a comparison of unimpacted animals for reference.

As with the slope plots, we need to prepare the data, and then plot it.

Data Preparation

To make this plot, we need estimates of health (anomaly) for each of the 6 injury classes for reproductively active and for non-reproductively active whales.

A reminder that this is the order of the categories:

gearInj | Definition -----------|---------- 1 | Severe with gear 2 | Moderate with gear 3 | Severe no gear 4 | Minor with gear 5 | Moderate no gear 6 | Minor no gear

We then take the tangRepro and tangNonRepro data frames that are created in the prepAmyEntanglementData.R script as part of tangled package construction. Within each of these data frames, we loop over the events and extract monthly estimates of the health anomaly during the event. Note that we also need a base reference class of whales that are not impacted by entanglement. We do that with the returnUnimpactedHealth.R function that is a helper function called from within prepBoxplotHealthData.R:

knitr::opts_chunk$set(warning=FALSE, message=FALSE)

library(tangled)
library(ggplot2)
library(dplyr)
library(magrittr)
library(stringr)

tmp <- prepBoxplotHealthData(tangRepro, tangNonRepro, anomFlag = FALSE, thold = 67)
dfLong <- rbind(tmp$dfLongRepro, tmp$dfLongNonRepro) 

At the end of this looping the new data look like this:

head(dfLong)

We can quickly summarise that to get a feel for how health changes as a function of gear injury type as well as an estimate of how many events are in each category. (Note that the number of events for the Unimpacted class does not reflect the actual number of events, rather the number of whale-months. So it's not useful to compare that to the entangled classes.)

dfLong %>% 
  group_by(status, gearInj) %>% 
  summarise(anomaly = median(hAnom, na.rm = TRUE))

dfLong %>% 
  filter(gearInj > 0) %>% 
  group_by(status, gearInj) %>% 
  summarise(numEvents = n())

Grouping by Decade

Goal here is to make the data groupings along the lines of the Meyer-Gutbrod paper that looks at the 80's, 90's, and 00's, i.e. prior to the big shift in distribution

dfLong <- dfLong %>% 
  mutate(decade = case_when(startDate >= 1980 & endDate < 1990  ~ "1980",
                            startDate >= 1990 & endDate < 2000  ~ "1990",
                            startDate >= 2000 & endDate < 2010  ~ "2000"),
         severity = case_when(gearInj == 6 ~ 'minor',
                              gearInj == 4 ~ 'minor',
                              gearInj == 5 ~ 'moderate',
                              gearInj == 2 ~ 'moderate',
                              gearInj == 3 ~ 'severe',
                              gearInj == 1 ~ 'severe',
                              gearInj == 0 ~ "Unimpacted")
         )

And then we want to summarize these by decade to get a feel for how much data we have in each group

dfLong %>% 
  filter(!is.na(decade)) %>% 
  filter(!is.na(gearInj)) %>% 
  group_by(decade, status, gearInj) %>% 
  summarise(numEvents = n())

ggplot(subset(dfLong, gearInj > 0 & !is.na(decade)))+
  geom_bar(aes(gearInj))+
  facet_grid(status ~ decade)

The Revised Plot

Ok, so then let's make the plot:

dfLong <- dfLong %>% 
  # filter(!is.na(decade)) %>% 
  filter(!is.na(gearInj))

pOut <- plotBoxplotHealth(dfLong, bsize = 12, cval = 3)
pOut
ggsave(plot = pOut, filename = 'boxplotHealth.pdf', path = '/Users/rob/Documents/research/manuscripts/KnowltonEtAl_Entanglement/images', device = 'pdf', width = 9, height = 6, units = 'in')
ggsave(plot = pOut, filename = 'boxplotHealth.png', path = '/Users/rob/Documents/research/manuscripts/KnowltonEtAl_Entanglement/images', device = 'png', dpi = 300, width = 9, height = 6, units = 'in')

Statistical Analysis of the Differences

With the data prepared and plotted, we want to determine quantitatively if there are differences between the entanglement classes and in between the reproductive status. To test this we can run a glm() and/or a glmm() on the dataset using health anomaly as the dependent variable.

library(multcomp)
dfglm <- transform(dfLong, gfac = as.factor(gearInj))
ft1 <- glm(hAnom ~ gfac, data = dfglm)
summary(ft1)
summary(glht(ft1, mcp(gfac = "Tukey")))

Include Reproductive Status

That shows some clear and significant patterns in the data. However, it doesn't show the comparison of the effects of reproductive status in the analysis. To do that, we'll explore two models: 1) a linear model with both injury and status as fixed effects; and 2) a mixed model with a fixed effect for injury, and then a random effect for reproductive status.

library(broom)
ft2 <- glm(hAnom ~ gfac + gfac:status, data = dfglm)
summary(ft2)
knitr::kable(tidy(ft2), digits = 3)

Ok, that was the glm and now we can look at the mixed effects model with a random effect

library(lme4)
library(sjPlot)
ft3 <- lmer(hAnom ~ gfac + (1 | status), data = dfglm)
summary(ft3)
sjp.lmer(ft3)

Summary Table

Need to put together a summary table for the ms of all the events included in this analysis.

df <- dfglm %>% 
  filter(gearInj != 0) %>% 
  group_by(gearInj, status) %>% 
  summarise(Total = n(),
            Length = mean(nMonths),
            maxLength = max(nMonths))
knitr::kable(df, digits = 2)

robschick/tangled documentation built on May 9, 2022, 4:07 p.m.