Introduction

Scenario: You have been asked to provide a preliminary stock assessment for Western Hake, a demersal species exploited by bottom trawlers on the Western Coast. The only available data is a 40-year catch time series (1984–2023). There is no survey index, no age data, and no length data.

From FishBase, you know the species has “Medium” resilience.

Your task: estimate MSY, evaluate current stock status (B/BMSY, F/FMSY), and provide management advice using the CMSY catch-only method (Froese et al. 2017).

Data file needed: catch_western_hake.csv

1 Setup

Install and load the required packages. The datalimited2 package implements CMSY and must be installed from GitHub.

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("datalimited2")) {
  if (!require("remotes")) install.packages("remotes")
  remotes::install_github("cfree14/datalimited2")
}

library(tidyverse)
library(datalimited2)

theme_set(theme_minimal(base_size = 14) +
            theme(panel.grid.minor = element_blank(),
                  plot.title = element_text(face = "bold")))

2 Load and Explore the Data

2.1 Read the data

catch_data <- read.csv("catch_western_hake.csv")
glimpse(catch_data)
#> Rows: 40
#> Columns: 2
#> $ year         <int> 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 199…
#> $ catch_tonnes <int> 19, 35, 46, 59, 71, 82, 88, 102, 105, 120, 130, 134, 143,…

2.2 Summary statistics

Before any modelling, we should understand what the data looks like:

cat("Time series:", min(catch_data$year), "–", max(catch_data$year),
    paste0("(", nrow(catch_data), " years)\n"))
#> Time series: 1984 – 2023 (40 years)
cat("Mean catch:  ", round(mean(catch_data$catch_tonnes)), "tonnes\n")
#> Mean catch:   69 tonnes
cat("Max catch:   ", round(max(catch_data$catch_tonnes)), "tonnes",
    paste0("(", catch_data$year[which.max(catch_data$catch_tonnes)], ")\n"))
#> Max catch:    143 tonnes (1996)
cat("Recent catch:", round(tail(catch_data$catch_tonnes, 1)), "tonnes",
    paste0("(", max(catch_data$year), ")\n"))
#> Recent catch: 26 tonnes (2023)

2.3 Plot the catch time series

This is all the data we have:

ggplot(catch_data, aes(year, catch_tonnes)) +
  geom_col(fill = "#1C7293", alpha = 0.8) +
  labs(title = "Western Hake — Catch Time Series",
       x = "Year", y = "Catch (tonnes)")
Western Hake catch time series. This is the only input to the CMSY model.

Western Hake catch time series. This is the only input to the CMSY model.

Q: Describe the catch trajectory in words. Can you identify distinct phases (development, peak, decline)?

Q: Does the shape of this catch history look “informative”? (Hint: a rise-then-fall pattern constrains what r and K can be. A flat catch would be consistent with many different stock dynamics.)

Q: Based on the catch pattern alone (no model), what would you guess about the current state of the stock?

3 Running CMSY

3.1 Background

The CMSY method (Froese et al. 2017) fits a Schaefer surplus production model using only catch data and a prior for the intrinsic growth rate r:

\[B_{t+1} = B_t + r \cdot B_t \left(1 - \frac{B_t}{K}\right) - C_t\]

It samples many (r, K) combinations, projects forward using observed catches, and filters out combinations that produce implausible biomass trajectories. The viable pairs define the posterior for reference points.

Resilience categories from FishBase map to r priors:

Resilience r prior range
High 0.6 – 1.5
Medium 0.2 – 0.8
Low 0.05 – 0.5
Very low 0.015 – 0.1

3.2 Fit the model

We run CMSY with “Medium” resilience — the appropriate category for Western Hake:

cmsy_fit <- cmsy2(
  year       = catch_data$year,
  catch      = catch_data$catch_tonnes,
  resilience = "Medium"
)
#> startbio= 0.2 0.6 default , intbio= 2019 0.01 0.4 default , endbio= 0.01 0.3 default 
#> First Monte Carlo filtering of r-k space with  10000  points...
#> 
#> Found  108  viable trajectories for 100  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.407 , 2.5 ]
#> Attempt  1  of  3  with  10000  additional points... 
#> 
#> Found altogether 396  viable trajectories for 345  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.312 , 2.47 ]
#> Attempt  2  of  3  with  20000  additional points... 
#> 
#> Found altogether 887  viable trajectories for 752  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.312 , 2.46 ]
#> Attempt  3  of  3  with  30000  additional points... 
#> 
#> Found altogether 1666  viable trajectories for 1377  r-k pairs
#> Final sampling in the tip area above r = 0.244 with 30000 additional points...
#> 
#> Found altogether 1922  viable trajectories for 1614  r-k pairs
#> ---------------------------------------
#> Catch data used from years 1984 - 2023 
#> Prior initial relative biomass = 0.2 - 0.6 default 
#> Prior intermediate rel. biomass= 0.01 - 0.4 in year 2019 default 
#> Prior final relative biomass   = 0.01 - 0.3 default 
#> Prior range for r = 0.2 - 0.8 default , prior range for k = 0.17 - 2.71 
#> 
#> Results of CMSY analysis 
#> -------------------------
#> Altogether 1922 viable trajectories for 1614  r-k pairs were found 
#> r   = 0.423 , 95% CL = 0.296 - 0.605 , k = 0.87 , 95% CL = 0.645 - 1.17 
#> MSY = 0.092 , 95% CL = 0.0816 - 0.104 
#> Relative biomass in last year = 0.117 k, 2.5th perc = 0.0149 , 97.5th perc = 0.291 
#> Exploitation F/(r/2) in last year = 1.2 
#> 
#> -------------------------------------------------------------
#> Fmsy = 0.212 , 95% CL = 0.148 - 0.303 (if B > 1/2 Bmsy then Fmsy = 0.5 r)
#> Fmsy = 0.0994 , 95% CL = 0.0695 - 0.142 (r and Fmsy are linearly reduced if B < 1/2 Bmsy)
#> MSY  = 0.092 , 95% CL = 0.0816 - 0.104 
#> Bmsy = 0.435 , 95% CL = 0.323 - 0.586 
#> Biomass in last year = 0.102 , 2.5th perc = 0.0129 , 97.5 perc = 0.253 
#> B/Bmsy in last year  = 0.235 , 2.5th perc = 0.0297 , 97.5 perc = 0.582 
#> Fishing mortality in last year = 0.255 , 2.5th perc = 0.103 , 97.5 perc = 2.01 
#> Exploitation F/Fmsy  = 2.56 , 2.5th perc = 1.03 , 97.5 perc = 20.2

3.3 Reference point estimates

ref_pts <- cmsy_fit$ref_pts
ref_ts  <- cmsy_fit$ref_ts

r_est   <- ref_pts$est[ref_pts$param == "r"]
K_est   <- ref_pts$est[ref_pts$param == "k"]
MSY_est <- ref_pts$est[ref_pts$param == "msy"]

cat("Estimated r   =", round(r_est, 3), "\n")
#> Estimated r   = 0.423
cat("Estimated K   =", round(K_est), "tonnes\n")
#> Estimated K   = 870 tonnes
cat("Estimated MSY =", round(MSY_est, 1), "tonnes\n")
#> Estimated MSY = 92 tonnes
ref_pts %>% knitr::kable(digits = 3, caption = "CMSY reference point estimates with 95% confidence intervals.")
CMSY reference point estimates with 95% confidence intervals.
param est lo hi
r 0.423 0.296 0.605
k 869.835 645.293 1172.510
msy 92.020 81.635 103.727
fmsy 0.212 0.148 0.303
bmsy 434.918 322.647 586.255

3.4 Diagnostic plots

The datalimited2 package provides a built-in multi-panel diagnostic:

plot_dlm(cmsy_fit)
CMSY diagnostic plot: catch, viable r-K pairs, B/Bmsy and F/Fmsy trajectories, and Kobe plot.

CMSY diagnostic plot: catch, viable r-K pairs, B/Bmsy and F/Fmsy trajectories, and Kobe plot.

Q: Look at the B/BMSY panel. When did the stock drop below BMSY?

Q: Look at the F/FMSY panel. When was overfishing most severe?

Q: In the Kobe plot, which quadrant is the stock currently in? (Red = overfished + overfishing; Green = healthy)

3.5 Detailed time series plots (nicer looking!)

3.5.1 B/BMSY trajectory

ggplot(ref_ts, aes(year)) +
  geom_ribbon(aes(ymin = bbmsy_lo, ymax = bbmsy_hi),
              fill = "#1C7293", alpha = 0.2) +
  geom_line(aes(y = bbmsy), colour = "#1C7293", linewidth = 1) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "#E07A5F") +
  annotate("text", x = min(ref_ts$year) + 1, y = 1.05,
           label = "Bmsy", colour = "#E07A5F") +
  labs(title = "Western Hake — Estimated B/Bmsy",
       x = "Year", y = expression(B/B[MSY]))
Estimated B/Bmsy with 95% confidence interval. The dashed line marks B = Bmsy.

Estimated B/Bmsy with 95% confidence interval. The dashed line marks B = Bmsy.

3.5.2 F/FMSY trajectory

ggplot(ref_ts, aes(year)) +
  geom_ribbon(aes(ymin = ffmsy_lo, ymax = ffmsy_hi),
              fill = "#1C7293", alpha = 0.2) +
  geom_line(aes(y = ffmsy), colour = "#1C7293", linewidth = 1) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "#E07A5F") +
  annotate("text", x = min(ref_ts$year) + 1, y = 1.05,
           label = "Fmsy", colour = "#E07A5F") +
  labs(title = "Western Hake — Estimated F/Fmsy",
       x = "Year", y = expression(F/F[MSY]))
Estimated F/Fmsy with 95% confidence interval. The dashed line marks F = Fmsy.

Estimated F/Fmsy with 95% confidence interval. The dashed line marks F = Fmsy.

4 Current Stock Status & Management Advice

final <- ref_ts %>% filter(year == max(year))

cat("Stock status in", max(ref_ts$year), ":\n\n")
#> Stock status in 2023 :
cat("  B/Bmsy =", round(final$bbmsy, 2),
    paste0(" [95% CI: ", round(final$bbmsy_lo, 2), " – ", round(final$bbmsy_hi, 2), "]\n"))
#>   B/Bmsy = 0.23  [95% CI: 0.03 – 0.58]
cat("  F/Fmsy =", round(final$ffmsy, 2),
    paste0(" [95% CI: ", round(final$ffmsy_lo, 2), " – ", round(final$ffmsy_hi, 2), "]\n"))
#>   F/Fmsy = 2.56  [95% CI: 1.03 – 20.25]
cat("  MSY    =", round(MSY_est, 1), "tonnes",
    paste0(" [95% CI: ", round(ref_pts$lo[ref_pts$param == "msy"], 1),
           " – ", round(ref_pts$hi[ref_pts$param == "msy"], 1), "]\n"))
#>   MSY    = 92 tonnes  [95% CI: 81.6 – 103.7]

4.0.1 Diagnosis (Stock status)

if (final$bbmsy < 1 & final$ffmsy > 1) {
  cat("→ OVERFISHED and OVERFISHING occurring (red Kobe quadrant)\n")
} else if (final$bbmsy < 1 & final$ffmsy <= 1) {
  cat("→ OVERFISHED but overfishing not currently occurring\n")
} else if (final$bbmsy >= 1 & final$ffmsy > 1) {
  cat("→ Not overfished but OVERFISHING occurring\n")
} else {
  cat("→ Stock in GOOD condition (green Kobe quadrant)\n")
}
#> → OVERFISHED and OVERFISHING occurring (red Kobe quadrant)

Q: Based on CMSY, what catch level would you recommend for next year? Should it be above, at, or below MSY? Why?

Q: How confident are you in this advice? (hint: Look at the width of the 95% CIs for MSY and B/BMSY.

5 Sensitivity Analysis

CMSY results depend heavily on the priors. A responsible assessment always tests: “How much do my results change if I change my assumptions?”

5.1 Sensitivity to the resilience prior

We re-run CMSY with different resilience assumptions:

cmsy_low    <- cmsy2(year = catch_data$year, catch = catch_data$catch_tonnes,
                      resilience = "Low")
#> startbio= 0.2 0.6 default , intbio= 2019 0.01 0.4 default , endbio= 0.01 0.3 default 
#> First Monte Carlo filtering of r-k space with  10000  points...
#> 
#> Found  1982  viable trajectories for 1287  r-k pairs
#> Final sampling in the tip area above r = 0.0826 with 10000 additional points...
#> 
#> Found altogether 3572  viable trajectories for 2486  r-k pairs
#> ---------------------------------------
#> Catch data used from years 1984 - 2023 
#> Prior initial relative biomass = 0.2 - 0.6 default 
#> Prior intermediate rel. biomass= 0.01 - 0.4 in year 2019 default 
#> Prior final relative biomass   = 0.01 - 0.3 default 
#> Prior range for r = 0.05 - 0.5 default , prior range for k = 0.271 - 10.9 
#> 
#> Results of CMSY analysis 
#> -------------------------
#> Altogether 3572 viable trajectories for 2486  r-k pairs were found 
#> r   = 0.219 , 95% CL = 0.123 - 0.388 , k = 1.58 , 95% CL = 0.827 - 3.04 
#> MSY = 0.0867 , 95% CL = 0.0602 - 0.125 
#> Relative biomass in last year = 0.117 k, 2.5th perc = 0.0143 , 97.5th perc = 0.289 
#> Exploitation F/(r/2) in last year = 1.28 
#> 
#> -------------------------------------------------------------
#> Fmsy = 0.109 , 95% CL = 0.0617 - 0.194 (if B > 1/2 Bmsy then Fmsy = 0.5 r)
#> Fmsy = 0.0512 , 95% CL = 0.0289 - 0.0907 (r and Fmsy are linearly reduced if B < 1/2 Bmsy)
#> MSY  = 0.0867 , 95% CL = 0.0602 - 0.125 
#> Bmsy = 0.792 , 95% CL = 0.413 - 1.52 
#> Biomass in last year = 0.185 , 2.5th perc = 0.0227 , 97.5 perc = 0.459 
#> B/Bmsy in last year  = 0.234 , 2.5th perc = 0.0287 , 97.5 perc = 0.579 
#> Fishing mortality in last year = 0.14 , 2.5th perc = 0.0567 , 97.5 perc = 1.14 
#> Exploitation F/Fmsy  = 2.74 , 2.5th perc = 1.11 , 97.5 perc = 22.4
cmsy_narrow <- cmsy2(year = catch_data$year, catch = catch_data$catch_tonnes,
                      r.low = 0.20, r.hi = 0.40)
#> startbio= 0.2 0.6 default , intbio= 2019 0.01 0.4 default , endbio= 0.01 0.3 default 
#> First Monte Carlo filtering of r-k space with  10000  points...
#> 
#> Found  237  viable trajectories for 222  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.544 , 2.59 ]
#> Attempt  1  of  3  with  10000  additional points... 
#> 
#> Found altogether 892  viable trajectories for 766  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.541 , 2.56 ]
#> Attempt  2  of  3  with  20000  additional points... 
#> 
#> Found altogether 2224  viable trajectories for 1869  r-k pairs
#> ---------------------------------------
#> Catch data used from years 1984 - 2023 
#> Prior initial relative biomass = 0.2 - 0.6 default 
#> Prior intermediate rel. biomass= 0.01 - 0.4 in year 2019 default 
#> Prior final relative biomass   = 0.01 - 0.3 default 
#> Prior range for r = 0.2 - 0.4 expert, , prior range for k = 0.339 - 2.71 
#> 
#> Results of CMSY analysis 
#> -------------------------
#> Altogether 2224 viable trajectories for 1869  r-k pairs were found 
#> r   = 0.331 , 95% CL = 0.278 - 0.394 , k = 1.07 , 95% CL = 0.845 - 1.36 
#> MSY = 0.0889 , 95% CL = 0.0759 - 0.104 
#> Relative biomass in last year = 0.115 k, 2.5th perc = 0.0159 , 97.5th perc = 0.29 
#> Exploitation F/(r/2) in last year = 1.27 
#> 
#> -------------------------------------------------------------
#> Fmsy = 0.166 , 95% CL = 0.139 - 0.197 (if B > 1/2 Bmsy then Fmsy = 0.5 r)
#> Fmsy = 0.0765 , 95% CL = 0.0642 - 0.091 (r and Fmsy are linearly reduced if B < 1/2 Bmsy)
#> MSY  = 0.0889 , 95% CL = 0.0759 - 0.104 
#> Bmsy = 0.537 , 95% CL = 0.422 - 0.682 
#> Biomass in last year = 0.124 , 2.5th perc = 0.0171 , 97.5 perc = 0.311 
#> B/Bmsy in last year  = 0.231 , 2.5th perc = 0.0318 , 97.5 perc = 0.58 
#> Fishing mortality in last year = 0.21 , 2.5th perc = 0.0835 , 97.5 perc = 1.52 
#> Exploitation F/Fmsy  = 2.74 , 2.5th perc = 1.09 , 97.5 perc = 19.9
get_summary <- function(fit, label) {
  rp <- fit$ref_pts
  ts <- fit$ref_ts %>% filter(year == max(year))
  tibble(
    Scenario   = label,
    r_est      = rp$est[rp$param == "r"],
    K_est      = rp$est[rp$param == "k"],
    MSY_est    = rp$est[rp$param == "msy"],
    MSY_lo     = rp$lo[rp$param == "msy"],
    MSY_hi     = rp$hi[rp$param == "msy"],
    BBmsy_now  = ts$bbmsy,
    FFmsy_now  = ts$ffmsy
  )
}

sensitivity <- bind_rows(
  get_summary(cmsy_fit,    "Medium (baseline)"),
  get_summary(cmsy_low,    "Low resilience"),
  get_summary(cmsy_narrow, "Manual r [0.2–0.4]")
)

sensitivity %>%
  mutate(across(where(is.numeric), ~round(., 2))) %>%
  knitr::kable(caption = "Sensitivity of CMSY results to the resilience prior.")
Sensitivity of CMSY results to the resilience prior.
Scenario r_est K_est MSY_est MSY_lo MSY_hi BBmsy_now FFmsy_now
Medium (baseline) 0.42 869.84 92.02 81.63 103.73 0.23 2.56
Low resilience 0.22 1584.92 86.66 60.24 124.67 0.23 2.74
Manual r [0.2–0.4] 0.33 1073.92 88.92 75.91 104.16 0.23 2.74 
ggplot(sensitivity, aes(Scenario, MSY_est)) +
  geom_col(fill = "#1C7293", alpha = 0.7, width = 0.6) +
  geom_errorbar(aes(ymin = MSY_lo, ymax = MSY_hi), width = 0.2) +
  labs(title = "MSY Sensitivity to Resilience Prior",
       x = "", y = "MSY (tonnes)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))
MSY estimates under different resilience priors. Error bars show 95% CIs.

MSY estimates under different resilience priors. Error bars show 95% CIs.

Q: Which prior gives the highest MSY? The lowest?

Q1: Does the management conclusion (overfished or not) change across scenarios? If all scenarios agree, does that increase your confidence?

5.2 Sensitivity to the final depletion prior

The depletion prior at the end of the time series anchors the biomass reconstruction. We now will test two extremes:

cmsy_depleted <- cmsy2(year = catch_data$year, catch = catch_data$catch_tonnes,
                        resilience = "Medium",
                        endb.low = 0.01, endb.hi = 0.20)
#> startbio= 0.2 0.6 default , intbio= 2019 0.01 0.4 default , endbio= 0.01 0.2 expert 
#> First Monte Carlo filtering of r-k space with  10000  points...
#> 
#> Found  69  viable trajectories for 66  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.334 , 2.59 ]
#> Attempt  1  of  3  with  10000  additional points... 
#> 
#> Found altogether 256  viable trajectories for 229  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.334 , 2.53 ]
#> Attempt  2  of  3  with  20000  additional points... 
#> 
#> Found altogether 623  viable trajectories for 558  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.334 , 2.5 ]
#> Attempt  3  of  3  with  30000  additional points... 
#> 
#> Found altogether 1201  viable trajectories for 1070  r-k pairs
#> Final sampling in the tip area above r = 0.239 with 30000 additional points...
#> 
#> Found altogether 1438  viable trajectories for 1300  r-k pairs
#> ---------------------------------------
#> Catch data used from years 1984 - 2023 
#> Prior initial relative biomass = 0.2 - 0.6 default 
#> Prior intermediate rel. biomass= 0.01 - 0.4 in year 2019 default 
#> Prior final relative biomass   = 0.01 - 0.2 expert 
#> Prior range for r = 0.2 - 0.8 default , prior range for k = 0.17 - 2.71 
#> 
#> Results of CMSY analysis 
#> -------------------------
#> Altogether 1438 viable trajectories for 1300  r-k pairs were found 
#> r   = 0.411 , 95% CL = 0.289 - 0.584 , k = 0.898 , 95% CL = 0.665 - 1.21 
#> MSY = 0.0922 , 95% CL = 0.08 - 0.106 
#> Relative biomass in last year = 0.0891 k, 2.5th perc = 0.0146 , 97.5th perc = 0.19 
#> Exploitation F/(r/2) in last year = 1.58 
#> 
#> -------------------------------------------------------------
#> Fmsy = 0.205 , 95% CL = 0.144 - 0.292 (if B > 1/2 Bmsy then Fmsy = 0.5 r)
#> Fmsy = 0.0732 , 95% CL = 0.0515 - 0.104 (r and Fmsy are linearly reduced if B < 1/2 Bmsy)
#> MSY  = 0.0922 , 95% CL = 0.08 - 0.106 
#> Bmsy = 0.449 , 95% CL = 0.332 - 0.607 
#> Biomass in last year = 0.0801 , 2.5th perc = 0.0131 , 97.5 perc = 0.17 
#> B/Bmsy in last year  = 0.178 , 2.5th perc = 0.0293 , 97.5 perc = 0.379 
#> Fishing mortality in last year = 0.325 , 2.5th perc = 0.153 , 97.5 perc = 1.98 
#> Exploitation F/Fmsy  = 4.44 , 2.5th perc = 2.09 , 97.5 perc = 27
cmsy_healthy <- cmsy2(year = catch_data$year, catch = catch_data$catch_tonnes,
                       resilience = "Medium",
                       endb.low = 0.40, endb.hi = 0.80)
#> startbio= 0.2 0.6 default , intbio= 2019 0.01 0.4 default , endbio= 0.4 0.8 expert 
#> First Monte Carlo filtering of r-k space with  10000  points...
#> 
#> Found  106  viable trajectories for 95  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.486 , 4.21 ]
#> Attempt  1  of  3  with  10000  additional points... 
#> 
#> Found altogether 383  viable trajectories for 319  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.464 , 3.18 ]
#> Attempt  2  of  3  with  20000  additional points... 
#> 
#> Found altogether 993  viable trajectories for 816  r-k pairs
#> Shrinking k space: repeating Monte Carlo in the interval [ 0.421 , 2.77 ]
#> Attempt  3  of  3  with  30000  additional points... 
#> 
#> Found altogether 1943  viable trajectories for 1582  r-k pairs
#> Final sampling in the tip area above r = 0.251 with 30000 additional points...
#> 
#> Found altogether 2238  viable trajectories for 1855  r-k pairs
#> ---------------------------------------
#> Catch data used from years 1984 - 2023 
#> Prior initial relative biomass = 0.2 - 0.6 default 
#> Prior intermediate rel. biomass= 0.01 - 0.4 in year 2019 default 
#> Prior final relative biomass   = 0.4 - 0.8 expert 
#> Prior range for r = 0.2 - 0.8 default , prior range for k = 0.339 - 8.14 
#> 
#> Results of CMSY analysis 
#> -------------------------
#> Altogether 2238 viable trajectories for 1855  r-k pairs were found 
#> r   = 0.445 , 95% CL = 0.278 - 0.711 , k = 0.836 , 95% CL = 0.611 - 1.14 
#> MSY = 0.0929 , 95% CL = 0.0824 - 0.105 
#> Relative biomass in last year = 0.533 k, 2.5th perc = 0.413 , 97.5th perc = 0.681 
#> Exploitation F/(r/2) in last year = 0.262 
#> 
#> -------------------------------------------------------------
#> Fmsy = 0.222 , 95% CL = 0.139 - 0.356 (if B > 1/2 Bmsy then Fmsy = 0.5 r)
#> Fmsy = 0.222 , 95% CL = 0.139 - 0.356 (r and Fmsy are linearly reduced if B < 1/2 Bmsy)
#> MSY  = 0.0929 , 95% CL = 0.0824 - 0.105 
#> Bmsy = 0.418 , 95% CL = 0.305 - 0.572 
#> Biomass in last year = 0.446 , 2.5th perc = 0.345 , 97.5 perc = 0.569 
#> B/Bmsy in last year  = 1.07 , 2.5th perc = 0.827 , 97.5 perc = 1.36 
#> Fishing mortality in last year = 0.0584 , 2.5th perc = 0.0457 , 97.5 perc = 0.0753 
#> Exploitation F/Fmsy  = 0.262 , 2.5th perc = 0.206 , 97.5 perc = 0.338
depl_sens <- bind_rows(
  get_summary(cmsy_fit,      "Default depletion"),
  get_summary(cmsy_depleted, "Assume very depleted"),
  get_summary(cmsy_healthy,  "Assume healthy")
)

depl_sens %>%
  select(Scenario, MSY_est, BBmsy_now, FFmsy_now) %>%
  mutate(across(where(is.numeric), ~round(., 2))) %>%
  knitr::kable(caption = "Sensitivity to assumed final depletion level.")
Sensitivity to assumed final depletion level.
Scenario MSY_est BBmsy_now FFmsy_now
Default depletion 92.0 0.23 2.56
Assume very depleted 92.2 0.18 4.44
Assume healthy 93.0 1.07 0.26

Q: How much does the final depletion assumption change the stock status conclusion? Is this reassuring or concerning?

Q: In a real assessment, what information would help you set the final depletion prior? (Think: expert knowledge, CPUE trends, survey snapshots)

Exercises

Exercise A — Management summary

Imagine you must brief a fisheries manager. Write a short paragraph (4–5 sentences) stating:

  1. The estimated MSY and its uncertainty
  2. Current stock status (B/BMSY, F/FMSY)
  3. Your catch recommendation for next year
  4. Your level of confidence in this advice (with justification)

Exercise B — The “declining catch” trap

Your colleague looks at the catch data and says: “Catches have dropped 40% since the peak — clearly the stock is collapsing.”

Do you agree? What alternative explanation exists? (Hint: catch might reflects BOTH stock abundance AND management decisions.)

Exercise C — Time series length

Run CMSY on different windows of the data:

# First 15 years only
cmsy_early <- cmsy2(year = catch_data$year[1:15],
                     catch = catch_data$catch_tonnes[1:15],
                     resilience = "Medium")

# Last 15 years only
cmsy_late <- cmsy2(year = catch_data$year[26:40],
                    catch = catch_data$catch_tonnes[26:40],
                    resilience = "Medium")

How do results compare to the full 40 years? Which time window is most informative and why?

Exercise D (advanced) — Compare catch-only methods

The datalimited2 package includes several other catch-only methods:

Method Function Reference
CMSY cmsy2() Froese et al. 2017
Bayesian surplus production bsm() Froese et al. 2017
Boosted regression tree zbrt() Zhou et al. 2017a
Optimized catch-only model ocom() Zhou et al. 2017b
Refined ORCS rorcs() Free et al. 2017

Try running one additional method and compare its B/BMSY estimate to CMSY. Do the methods agree?