Summary

To calculate the status of global fish stocks from the SAUP 2022 catch data we use catch-MSY (CMSY) developed by Martell and Froese (2012). This model requires catch and resilience information for each stock. Stocks are defined by FAO areas and are limited to only those records where catch is reported at the species level.

The datalimited R package, developed by Sean Anderson, is used to run CMSY.

Updates from previous scenario year

None. This script wasn’t run in v2023. Associated data files within the repo were copied from v2022 to v2023.

Setup

To install the datalimited package:

install.packages("remotes")
remotes::install_github("datalimited/datalimited")

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',message = FALSE, warning = FALSE)

#devtools::install_github("datalimited/datalimited")
library(datalimited) #has the 4 catch only models
library(tidyverse)
library(doParallel)
library(here)

setwd(here::here("globalprep/fis/v2022"))


source(here('workflow/R/common.R'))

Load catch data

Read in catch data aggregated from half degree cell to OHI region catch_data_prep.Rmd.

catch <- read_csv('output/stock_catch.csv') %>%
          rename(common = CommonName)

Catch-MSY

If the CMSY values have already been calculated for each stock, this step can be skipped:

registerDoParallel(cores = 5)

# calculate cmsy values and create a .rds file that is saved on Mazu

cmsy_fits <- plyr::dlply(catch, c("stock_id", "common"), function(x) {
  
    #make sure the data is ordered from 1950 to 2019
    x <- arrange(x,year)
    
    out <- cmsy(ct = x$tons, yr = x$year,  start_r = resilience(x$Resilience[1]), 
      reps = 2e4)
    
    out$year <- x$year
    out
    
  }, .parallel = TRUE)


saveRDS(cmsy_fits, file = file.path(dir_M,"git-annex/globalprep/fis/v2022/int/cmsy-fits.rds"))

Take output and format:

cmsy_fits <- readRDS(file.path(dir_M,"git-annex/globalprep/fis/v2022/int/cmsy-fits.rds"))

fake_data <- data.frame(bbmsy_q2.5 = NA, bbmsy_q25 = NA, bbmsy_q50 = NA, 
  bbmsy_q75 = NA, bbmsy_q97.5 = NA)

 cmsy_bbmsy <- plyr::ldply(cmsy_fits, function(x) {
  bbmsy_cmsy <- x$biomass[, -1] / x$bmsy
  bbmsy_out <- tryCatch({
    bbmsy_out <- summarize_bbmsy(bbmsy_cmsy)
    bbmsy_out$year <- x$year
    bbmsy_out}, error = function(e) fake_data)
})

cmsy_bbmsy$model <- "CMSY"

write.csv(cmsy_bbmsy,file='output/cmsy_bbmsy.csv', row.names=FALSE)

Explore why there are some NAs (I think non convergance)

cmsy_bbmsy <- read_csv("output/cmsy_bbmsy.csv") 

nas <- cmsy_bbmsy %>%
  group_by(stock_id)%>%
  summarize(m = mean(bbmsy_mean))%>%
  filter(is.na(m))

nrow(nas)

Results

Looking at mean bbmsy for 2014

results <- read_csv('output/cmsy_bbmsy.csv') %>%
            filter(year == 2014)
hist(results$bbmsy_mean)

Why are so many around 0.5 (extracting half of maximum sustainable yield)?

stocks_0.5 <- results %>%
              filter(bbmsy_mean < 0.6,
                     bbmsy_mean > 0.4)

Anderson, S.C., Cooper, A.B., Jensen, O.P., et al. (2017) Improving estimates of population status and trend with superensemble models. Fish and Fisheries.

OHI 2023: Food Provision, Calculating fisheries stock status