WGMMUME Climate Change subgroup - California sea lion example

Author

Len Thomas

Published

February 5, 2025

library(readxl)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ purrr::%||%()   masks base::%||%()
✖ 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

Introduction

Here we present a preliminary analysis of California sea lion (CSL) strandings data provided by Justin Greenman on 22nd November 2024.

We also make use of population estimates provided by Laake et al. (2008), which (so far as I currently know) are the most recent for CSL.

Read in the data

#read in estimated population sizes from Laake paper
pop_size <- read.delim("Laake_2018.txt", sep = " ")
#Just keep 1st and 5th columns
pop_size <- pop_size[, c(1, 5)]
colnames(pop_size) <- c("year", "N")
#Strip of "a" from some years and make numeric
pop_size$year <- as.numeric(substr(pop_size$year, 1, 4))
#Strip any comma from N and make numeric
pop_size$N <- as.numeric(gsub(",", "", pop_size$N))
plot(pop_size$year, pop_size$N, type = "l", xlab = "year", ylab = "N")

#read in strandings data
infile <- "WCR Zc-At 1982-2023.xlsx"
tab1 <- read_xlsx(infile, sheet = 1)
tab1$year <- year(tab1$"Date of Occurance")
tab1$dead <- tab1$"Initial Condition" %in% c("FRESH DEAD", "MODERATE DECOMPOSITION",
                                       "ADVANCED DECOMPOSITION")
#Could get further information by looking for "DIED" in "Final Disposition remarks" and "Stranding Remarks" columns
tab2 <- read_xlsx(infile, sheet = 2)
tab2$year <- tab2$"Year of Observation"
#dead excludes "unknown" as well as "alive"
tab2$dead <- tab2$"Observation Status" %in% c("Fresh Dead", "Moderate Decomposition",
                                               "Advanced Decomposition", "Mummified/Skeletal")

tab3 <- read_xlsx(infile, sheet = 3)
tab3$year <- tab3$"Year of Observation"
tab3$dead <- tab3$"Observation Status" %in% c("Fresh Dead", "Moderate Decomposition",
                                              "Advanced Decomposition", "Mummified/Skeletal")

all_strandings <- rbind(tab1[, c("year", "dead")], tab2[, c("year", "dead")], tab3[, c("year", "dead")])
strandings <- all_strandings %>% group_by(year) %>%
  summarize(stranded = n(), dead = sum(dead)) %>%
  full_join(pop_size[, c("year", "N")], by = "year") %>%
  mutate(p_stranded = stranded / N, p_dead = dead / N) %>%
  arrange(year)

Data summary

Below is a graphical summary of the strandings data.

par(mfrow = c(2, 2))
#with(strandings, plot(year, N, type = "l", main = "Estimated pop size"))
#plot(1:10, 1:10, axes = FALSE, type = "n", xlab = "", ylab = "")
with(strandings, plot(year, stranded, type = "l", main = "Recorded strandings"))
with(strandings, plot(year, dead, type = "l", main = "Recorded deaths"))
with(strandings, plot(year, p_stranded, type = "l", main = "Proportion of pop stranded"))
with(strandings, plot(year, p_dead, type = "l", main = "Proportion of pop dead"))

par(mfrow = c(1, 1))

Recorded strandings and deaths have been going up, but when you divide by estimated population size, there is a flat trend over time. (Note that the time span of the strandings data and population estimates doesn’t fully overlap, so the bottom plots are for fewer time steps.)

UME detection

In the plot below, I show recorded dead stranded animals as a black line, the mean of the previous 5 year’s deaths as the red line, +/- 2 SD as red dashed lines, and strandings that are outside the red dashed lines as red dots.

#Work out the mean of the last 5 years, and +/- 2SD
strandings$mv_mean_dead <- strandings$mv_mean_pdead <- strandings$sd_mv_mean_dead <- 
  strandings$lcl_mv_mean_dead <- strandings$ucl_mv_mean_dead <- NA
lag <- 5
for(i in (lag + 1):nrow(strandings)){
  strandings$mv_mean_dead[i] <- mean(strandings$dead[(i-(lag+1)):(i-1)])
  strandings$sd_mv_mean_dead[i] <- sd(strandings$dead[(i-(lag+1)):(i-1)])
  strandings$lcl_mv_mean_dead[i] <- strandings$mv_mean_dead[i] - 2 * strandings$sd_mv_mean_dead[i] 
  strandings$ucl_mv_mean_dead[i] <- strandings$mv_mean_dead[i] + 2 * strandings$sd_mv_mean_dead[i] 
}

with(strandings, plot(year, dead, type = "l", main = "Recorded deaths"))
with(strandings, lines(year, mv_mean_dead, col = "red"))
with(strandings, lines(year, lcl_mv_mean_dead, col = "red", lty = 2))
with(strandings, lines(year, ucl_mv_mean_dead, col = "red", lty = 2))
ind <- (strandings$dead > strandings$ucl_mv_mean_dead) | (strandings$dead < strandings$lcl_mv_mean_dead)
ind[is.na(ind)] <- FALSE
with(strandings, points(year[ind], dead[ind], col = "red", pch = 19))

We were outside the +/- 2SD bounds in the following 6 years:

strandings$year[ind]

[1] 1991 1992 1998 2009 2022 2023

A real-world algorithm would likely discount previous unusual values in calculating the mean, and we may use more than 5 years.

Discussion

In general, it appears the above approach at least partly deals with the fact that the recorded deaths are trending upwards (in this case because the population is trending upwards). So, any kind of gradual increase in deaths over time will end not to be declared a UME.

Some refinements would be to fit a smooth to population size (and perhaps project forward a few years), to discount previous unusual values (see previous section), to adjust the lag period, to work on the log scale so that we don’t get negative lower bounds … there are likely lots of other possible refinements!

References

Laake, J.L., M.S. Lowry, R.L. DeLong, S.R. Melin and J.V. Carretta. 2018. Population growth and status of California sea lions. Journal of Wildlife Management 82: 583-595.