Skip to contents

Note: The model-fitting chunks in this vignette are disabled because of a Matrix compatibility issue. See the README for details.

library(magrittr) # required for this vignette to work

1. Install lbmstoolbox

To run this vignette, install lbmstoolbox in one of R’s library paths. Install it system-wide with devtools::install_github() or locally with renv::install(). When running the source interactively, use the vignettes directory as the working directory so that the bundled data file can be found.

2. Read sampling and biological data

See the LB-SPR tutorial for details of the required input data.

spp_sampling_details <- readRDS("data/spp_test_data.rds")
catch_data <- spp_sampling_details$catch_data

3. Run LIME model

3.1 Build biological parameters for LIME

See the LIME documentation for parameter definitions.

The example biological-parameter list contains fields for both LB-SPR and LIME. Remove the LB-SPR-specific field. rec_variability_mean and rec_variability_sd define the mean and standard deviation of LIME’s recruitment-variability prior (SigRprior), preventing unrealistically variable recruitment estimates (Rudd and Thorson, 2017). The example values are based on Thorson et al. (2014).

bio_params <- spp_sampling_details$bio_params
bio_params$fecb <- NULL
bio_params
#> $linf
#> [1] 42.41
#> 
#> $k
#> [1] 0.21
#> 
#> $t0
#> [1] -1.202
#> 
#> $lwa
#> [1] 0.012
#> 
#> $lwb
#> [1] 3.011
#> 
#> $l50
#> [1] 19.45
#> 
#> $M
#> [1] 0.354
#> 
#> $max_age
#> [1] 9
#> 
#> $rec_variability_mean
#> [1] 0.748
#> 
#> $rec_variability_sd
#> [1] 0.293
#> 
#> $max_age
#> [1] 9
#> 
#> $l95
#> [1] 22.37

3.2 Build exploitation parameters for LIME

Unlike LB-SPR, LIME requires selectivity priors and a penalty that constrains year-to-year changes in estimated fishing mortality (Rudd and Thorson, 2017). Here, SL50 is halfway between the smallest observed length class and the modal class. SL95 is set to SL50 * 1.15, following the approach used by Prince et al. (2022) for the maturation ogive. sigmaF uses the default value of 0.2 from Rudd and Thorson (2017).

summary_catch <- catch_data$long %>%
  dplyr::select(MeanLength, catch) %>%
  dplyr::group_by(MeanLength) %>%
  dplyr::summarise(total_catch = sum(catch))
catch_mode <- summary_catch$MeanLength[which.max(summary_catch$total_catch)]
sl50 <- median(min(summary_catch$MeanLength):catch_mode)
sel_params <- list(
  sl50 = sl50,
  sl95 = sl50 * 1.15
)

exploitation_params <- list(
  s50 = sel_params$sl50,
  s95 = sel_params$sl95,
  sigmaF = 0.2
)
exploitation_params
#> $s50
#> [1] 12.5
#> 
#> $s95
#> [1] 14.375
#> 
#> $sigmaF
#> [1] 0.2

3.3 Build the remaining LIME context

As an age-structured model, LIME must be told whether selectivity and maturity are supplied by length or age. This list also sets steepness to h = 0.7.

extra_lht_details <- list(
  binwidth = 1,
  selex_input = "length",
  maturity_input = "length",
  nseasons = 1,
  nfleets = 1,
  h = 0.7
)

3.4 Instantiate and run LIME

With simul = TRUE, the wrapper returns the following additional output for each time step:

  • Observed catch-at-length
  • Observed catch-at-length, scaled to sum to one
  • Standardised fished and unfished population-at-age structures
  • Standardised estimated catch-at-length after applying estimated selectivity to the fished population
  • Standardised estimated catch-at-length after applying estimated selectivity to the unfished population
# Find SL50 and SL95 prior
iim_lime_logistic <- lbmstoolbox::LimeLbms$new(
  bio_params,
  exploitation_params,
  catch_data
)
lht_context <- iim_lime_logistic$build_lht_context(extra_lht_details)
result_lime <- iim_lime_logistic$run(lht_context, simul = TRUE)

3.5 Breaking down results

Results have two components: evaluation, which contains model estimates, and simulation, which contains the additional output. The latter also retains the raw LIME output.

The estimates data frame also includes confidence intervals for F, SPR, and recruitment (R), with lower_ and upper_ prefixes.

3.5.1 Estimates dataframe

result_lime$evaluation$estimates %>%
  head(3)

3.5.2 Simulation dataframes

The simulation output contains separate catch and population data frames. Their relevant columns are:

  • unfished_pop: age structure of the unfished population
  • fished_pop: age structure of the fished population
  • unfished_vul: estimated catch from the unfished population after applying selectivity
  • exp_catch: estimated catch from the fished population after applying selectivity
  • real_catch: observed catch in absolute numbers
  • relative_catch: observed catch scaled to sum to one

Unlike LB-SPR, LIME returns catch and population structures separately because catch is described by length and population by age.

result_lime$simulation$model_info$catch %>%
  head(3)
result_lime$simulation$model_info$population %>%
  head(3)

4. Plot examples

4.1 Estimates plot

A quick example on how to plot SPR and its confidence interval is shown below.

estimates <- result_lime$evaluation$estimates
g <- estimates %>%
  ggplot2::ggplot(ggplot2::aes(x = years, y = SPR)) +
  ggplot2::geom_point() +
  ggplot2::geom_line() +
  ggplot2::geom_ribbon(ggplot2::aes(x = years, ymin = lower_spr, ymax = upper_spr, fill = "maroon2", alpha = 0.1)) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "none")

g

4.2 Simulation plot

The first plot shows estimated catch-at-length for the fished population in 2015. The second compares the fished and unfished population-at-age structures for the same year.

4.2.1 Estimated catch

catch_2015 <- result_lime$simulation$model_info$catch %>%
  dplyr::filter(year == 2015) %>%
  dplyr::select(lengths, exp_catch)

g <- catch_2015 %>%
  ggplot2::ggplot(ggplot2::aes(x = lengths, y = exp_catch)) +
  ggplot2::geom_point() +
  ggplot2::geom_line() +
  ggplot2::xlab("Lengths") +
  ggplot2::ylab("Standardised relative catch")

g

4.2.2 Estimated populations

population_2015 <- result_lime$simulation$model_info$population %>%
  dplyr::filter(year == 2015) %>%
  dplyr::select(age, unfished_pop, fished_pop) %>%
  tidyr::pivot_longer(!age, names_to = "variables", values_to = "values")

g <- population_2015 %>%
  ggplot2::ggplot(ggplot2::aes(x = age, y = values, color = variables)) +
  ggplot2::geom_point() +
  ggplot2::geom_line() +
  ggplot2::xlab("Ages") +
  ggplot2::ylab("Standardised relative populations") +
  ggplot2::scale_color_manual(
    name = "Catch Type",
    values = c(
      unfished_pop = "seagreen",
      fished_pop = "tomato"
    )
  )

g

5. References

  1. Rudd, M.B., Thorson, J.T., 2017. Accounting for Variability and Biases in Data-limited Fisheries Stock Assessment. Canadian Journal of Fisheries and Aquatic Sciences 75, 1019–1035. https://doi.org/10.1139/cjfas-2017-0143
  2. Thorson, J.T., Jensen, O.P., Zipkin, E.F., 2014. How variable is recruitment for exploited marine fishes? A hierarchical model for testing life history theory. Can. J. Fish. Aquat. Sci. 71, 973–983. https://doi.org/10.1139/cjfas-2013-0645
  3. Prince, J., Harford, W.J., Taylor, B.M., Lindfield, S.J., 2022. Standard histological techniques systematically under‐estimate the size fish start spawning. Fish and Fisheries 23, 1507–1516. https://doi.org/10.1111/faf.12702