Working with Recruitment Deviates • marlin

Fish populations are often characterized by variation, particularly in the process known as “recruitment”, the entry of new individuals to the population through reproduction (as opposed to immigration).

marlin allows users a range of options for simulating variability in the recruitment process. At the most basic, users can specify the variation of “recruitment deviations”, deviations in the value of recruitment away from the mean level of recruitment, for example the number of recruits given spawning stock biomass (SSB) given a Beverton-Holt spawner-recruit relationship. From there, users can specify whether recruitment deviations are autocorrelated, and lastly the degree of correlation in recruitment deviates among species.

Users can also supply their own external recruitment deviates, allowing users to for example generate recruitment deviates that are correlated across space, time, and species in complex ways.

To start with, we’ll simulate two uncorrelated species with varying degrees of recruitment variation (sigma_rec) and autocorrelation in recruitment variation (ac_rec)

library(marlin)

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() ──
#> ✖ 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

library(marlin)

theme_set(marlin::theme_marlin())

resolution <- c(10,10) # resolution is in squared patches, so 20 implies a 20X20 system, i.e. 400 patches 

patch_area <- 10

seasons <- 2

years <- 20

tune_type <- "depletion"

steps <- years * seasons

yft_diffusion <- 2
  
yft_depletion <- .2

rec_factor <- 10

mako_depletion <- 0.3

mako_diffusion <- 100

critters <- c("yellowfin", "mako")



fauna <- 
  list(
    "yellowfin" = create_critter(
      scientific_name = "Thunnus albacares",
      adult_diffusion = yft_diffusion, # cells per year
      recruit_diffusion = rec_factor * yft_diffusion,
      fished_depletion = yft_depletion, # desired equilibrium depletion with fishing (1 = unfished, 0 = extinct),
      density_dependence = "pre_dispersal", # recruitment form, where 1 implies local recruitment
      seasons = seasons,
      ssb0 = 5000,
      sigma_rec = 0.42,
      ac_rec = 0.42
      ),
    "mako" = create_critter(
      scientific_name = "Isurus oxyrinchus",
      adult_diffusion = mako_diffusion,
      recruit_diffusion = rec_factor * .1,
      fished_depletion = mako_depletion,
      density_dependence = "local_habitat", # recruitment form, where 1 implies local recruitment
      burn_years = 200,
      ssb0 = 100,
      seasons = seasons,
      fec_form = "pups",
      pups = 6,
      sigma_rec = 0.01,
      ac_rec = 0
    )
  )

fleets <- list("longline" = create_fleet(list(
  `yellowfin` = Metier$new(
    critter = fauna$`yellowfin`,
    price = 100, # price per unit weight
    sel_form = "logistic", # selectivity form, one of logistic or dome
    sel_start = .3, # percentage of length at maturity that selectivity starts
    sel_delta = .1, # additional percentage of sel_start where selectivity asymptotes
    catchability = .01, # overwritten by tune_fleet but can be set manually here
    p_explt = 1
    ),
  `mako` = Metier$new(
    critter = fauna$`mako`,
    price = 10,
    sel_form = "logistic",
    sel_start = .1,
    sel_delta = .01,
    catchability = 0,
    p_explt = 1
  )),
  mpa_response = "stay",
  base_effort = 1000 * prod(resolution),
  spatial_allocation = "rpue",
  resolution = resolution
))

a <- Sys.time()

fleets <- tune_fleets(fauna, fleets, tune_type = tune_type) # tunes the catchability by fleet to achieve target depletion

Sys.time() - a
#> Time difference of 14.50491 secs

# run simulations

a <- Sys.time()

recruitment_sim <- simmar(fauna = fauna,
                  fleets = fleets,
                  steps = steps)

Sys.time() - a
#> Time difference of 0.05609894 secs

sim <- process_marlin(recruitment_sim)


sim$fauna |> 
  filter(age == min(age)) |> 
  ggplot(aes(step, n)) + 
  geom_line() + 
  geom_point() + 
  facet_wrap(~critter, scales = "free_y")

Timeline of recruits (age 0 fish) into the population over time given degrees of variation and autocorrelation in the recruitment process.

We can expand on this simple analysis by adding in a new feature, a correlation matrix for the recruitment deviates between the two simulated species. As this is a simple two-species system, this is just a two x two matrix with the off-diagonal elements indicating the correlation between the recruitment deviates of the two species. We will then pass this matrix to the simmar function.


critter_correlations <- matrix(c(1, -.8, -.8, 1), nrow = 2, byrow = TRUE)


fauna <- 
  list(
    "yellowfin" = create_critter(
      scientific_name = "Thunnus albacares",
      adult_diffusion = yft_diffusion, # cells per year
      recruit_diffusion = rec_factor * yft_diffusion,
      fished_depletion = yft_depletion, # desired equilibrium depletion with fishing (1 = unfished, 0 = extinct),
      density_dependence = "pre_dispersal", # recruitment form, where 1 implies local recruitment
      seasons = seasons,
      ssb0 = 5000,
      sigma_rec = 0.42,
      ac_rec = 0
      ),
    "mako" = create_critter(
      scientific_name = "Isurus oxyrinchus",
      adult_diffusion = mako_diffusion,
      recruit_diffusion = rec_factor * .1,
      fished_depletion = mako_depletion,
      density_dependence = "local_habitat", # recruitment form, where 1 implies local recruitment
      burn_years = 200,
      ssb0 = 100,
      seasons = seasons,
      fec_form = "pups",
      pups = 6,
      sigma_rec = 0.42,
      ac_rec = 0
    )
  )

fleets <- list("longline" = create_fleet(list(
  `yellowfin` = Metier$new(
    critter = fauna$`yellowfin`,
    price = 100, # price per unit weight
    sel_form = "logistic", # selectivity form, one of logistic or dome
    sel_start = .3, # percentage of length at maturity that selectivity starts
    sel_delta = .1, # additional percentage of sel_start where selectivity asymptotes
    catchability = .01, # overwritten by tune_fleet but can be set manually here
    p_explt = 1
    ),
  `mako` = Metier$new(
    critter = fauna$`mako`,
    price = 10,
    sel_form = "logistic",
    sel_start = .1,
    sel_delta = .01,
    catchability = 0,
    p_explt = 1
  )),
  mpa_response = "stay",
  base_effort = 1000 * prod(resolution),
  spatial_allocation = "rpue",
  resolution = resolution
))

a <- Sys.time()

fleets <- tune_fleets(fauna, fleets, tune_type = tune_type) # tunes the catchability by fleet to achieve target depletion

Sys.time() - a
#> Time difference of 6.139112 secs

# run simulations

a <- Sys.time()

recruitment_sim <- simmar(
  fauna = fauna,
  fleets = fleets,
  steps = steps,
  cor_rec = critter_correlations
)

Sys.time() - a
#> Time difference of 0.05300307 secs

sim <- process_marlin(recruitment_sim)


sim$fauna |> 
  filter(age == min(age)) |> 
  select(step,patch,critter, n) |> 
  pivot_wider(names_from = critter, values_from = n) |> 
  ggplot(aes(mako, yellowfin)) + 
  geom_point()

Scatter plot of negatively correlated recruitment between the two simulated species.

Users can also supply their own matrix of recruitment deviates. This matrix has number of columns equal to the number of critters, and rows equal to the number of years to be simulated plus one additional number of seasons per year (a required buffer). This step can be useful if for example the user wants to share the same vector of recruitment deviates across multiple simulations, or generate more complex recruitment deviates driven by for example environmental covariates.



  sigma_recs <- purrr::map_dbl(fauna, "sigma_rec") # gather recruitment standard deviations

  ac_recs <- purrr::map_dbl(fauna, "ac_rec") # gather autocorrelation in recruitment standard deviations

  n_critters <-  length(fauna)
  
  covariance_rec <- critter_correlations * (sigma_recs %o% sigma_recs)

  rec_steps <- steps + seasons
  
  log_rec_devs <- matrix(NA, nrow = rec_steps, ncol = n_critters, dimnames = list(1:(rec_steps), names(fauna)))

  log_rec_devs[1,] <- mvtnorm::rmvnorm(1,rep(0, n_critters),sigma = covariance_rec)

  for (i in 2:rec_steps){

    log_rec_devs[i, ] <- ac_recs *  log_rec_devs[i - 1, ] + sqrt(1 - ac_recs ^ 2) * mvtnorm::rmvnorm(1, rep(0,n_critters), sigma = covariance_rec)

  }
  
  a <- Sys.time()

recruitment_sim <- simmar(fauna = fauna,
                  fleets = fleets,
                  steps = steps,
                  log_rec_devs = log_rec_devs)

Sys.time() - a
#> Time difference of 0.1385801 secs

sim <- process_marlin(recruitment_sim)


sim$fauna |> 
  filter(age == min(age)) |> 
  select(step,patch,critter, n) |> 
  pivot_wider(names_from = critter, values_from = n) |> 
  ggplot(aes(mako, yellowfin)) + 
  geom_point()

Scatter plot of externally supplied negatively correlated recruitment between the two simulated species.