inst/markdown/archive/04.snowcrab_fishery_model_turing.md

In jae0/bio.snowcrab: Snow crab assessment suite using aegis* and associated projects and data streams.

title: "Logistic Fishery Model" author: "Jae S. Choi" toc: true number-sections: true highlight-style: pygments editor: render-on-save: false markdown: references: location: document format: html: code-fold: true html-math-method: katex self-contained: true embed-resources: true pdf: pdf-engine: lualatex docx: default

Purpose

Prepare aggregated time-series data for inference using a Bayesian fishery model. Currently, the default model is the discrete logistic formulation. It is labelled, "logistic_historical" in the following. There are other variations possible. See the calling code.

Prepare data for discrete logistic models

This step requires bio.snowcrab/inst/scripts/03.biomass_index_carstm.r to be completed. That step statistically models the survey samples in space and time and then aggregates across space for an estimate/index of abundance.

Here we read in these results and prepare (format) them for timeseries analysis.

#| eval: true
#| output: false
  require(aegis)
  year.assessment = 2023
  p = bio.snowcrab::load.environment( year.assessment=year.assessment )
  modeldir = p$modeldir
  snowcrab_filter_class = "fb"     # fishable biomass (including soft-shelled )  "m.mat" "f.mat" "imm"

  carstm_model_label = paste( "default", snowcrab_filter_class, sep="_" )
  fishery_model_label = "turing1"

  # save to : carstm_directory = file.path(modeldir, carstm_model_label)

  fishery_model_data_inputs( 
    year.assessment=year.assessment,   
    type="biomass_dynamics", 
    snowcrab_filter_class="fb",
    modeldir=modeldir, 
    carstm_model_label=carstm_model_label,
    for_julia=TRUE,
    fishery_model_label="turing1"  
  ) 

  # the name/location of the above creates a data file (which will be passed to Julia, later) 
  fndat  = file.path( modeldir, carstm_model_label, "biodyn_biomass.RData" )  # 

  # for development using dynamical_model/snowcrab, save/copy at alternate locations
  # modeldir = file.path( homedir, "projects", "dynamical_model", "snowcrab", "data" ) 

Model fitting and parameter inference

Now that data inputs are ready, we can continue on to model fitting and inference. This is done in (Julia)[https://julialang.org/].

First install or download Julia and install appropriate libraries from withing Julia.

There are two ways you can use this.

1. Run within Julia directly.

   This method is most flexible and documented here, where more options are also shown. But it will require learning the language Julia and Turing library, but is very simple.

#| eval: false 
#| output: false

  year_assessment =2023
  yrs = 1999:year_assessment

  rootdir = homedir()

  project_directory  = joinpath( rootdir, "projects", "bio.snowcrab", "inst", "julia" ) 
  bio_data_directory = joinpath( rootdir, "bio.data", "bio.snowcrab", "modelled", "default_fb" )
  outputs_directory  = joinpath( rootdir, "bio.data", "bio.snowcrab", "fishery_model" ) 

  model_variation = "logistic_discrete_historical"   
  model_outdir = joinpath( outputs_directory, string(year_assessment), model_variation )

  include( joinpath(project_directory, "startup.jl" ) )  # load libs
  include( joinpath(project_directory, "logistic_discrete_functions.jl" ) )  # load core modelling functions

  o = load( joinpath( bio_data_directory, "biodyn_biomass.RData" ), convert=true)   # load data; alter file path as required   

  # run the whole script or in parts inside "bio.snowcrab/inst/julia/logistic_discrete.jl" for more control
  aulab ="cfanorth"     
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

  aulab ="cfasouth"    
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

  aulab ="cfa4x"      
  include( joinpath(project_directory, "logistic_discrete.jl" ) )   

And that is it. Now we continue to look at the results:

Stock status: Biomass Index (aggregate) ... {.c}

\begin{tiny} ```{r fbindex-timeseries, out.width='60%', echo=FALSE, fig.align='center', fig.cap = 'The fishable biomass index (t) predicted by CARSTM of Snow Crab survey densities. Error bars represent Bayesian 95\% Credible Intervals. Note large errors in 2020 when there was no survey.' } include_graphics( file.path( SCD, 'modelled', 'default_fb', 'aggregated_biomass_timeseries' , 'biomass_M0.png') )

\@ref(fig:fbindex-timeseries)

\end{tiny}




## Stock status: Modelled Biomass (pre-fishery) ... {.c}

\begin{tiny}
```{r logisticPredictions, out.width='32%', echo=FALSE, fig.show='hold', fig.align='center', fig.cap = 'Model 1 fishable, posterior mean modelled biomass (pre-fishery; kt) are shown in dark orange for N-ENS, S-ENS and 4X (left, middle and right). Light orange are posterior samples of modelled biomass (pre-fishery; kt) to illustrate the variability of the predictions. The biomass index (post-fishery, except prior to 2004) after model adjustment by the model catchability coefficient is in gray.' } 
loc = file.path( SCD, 'fishery_model', year.assessment, 'logistic_discrete_historical' )
fn1 = file.path( loc, 'plot_predictions_cfanorth.pdf' ) 
fn2 = file.path( loc, 'plot_predictions_cfasouth.pdf' ) 
fn3 = file.path( loc, 'plot_predictions_cfa4x.pdf' ) 
include_graphics(c(fn1, fn2, fn3) )
# \@ref(fig:logisticPredictions)

\end{tiny}

Stock status: Fishing Mortality ... {.c}

```{r logisticFishingMortality, out.width='32%', echo=FALSE, fig.show='hold', fig.align='center', fig.cap = 'Time-series of modelled instantaneous fishing mortality from Model 1, for N-ENS (left), S-ENS (middle), and 4X (right). Samples of the posterior densities are presented, with the darkest line being the mean.' } odir = file.path( fishery_model_results, year.assessment, "logistic_discrete_historical" ) fn1 = file.path( odir, "plot_fishing_mortality_cfanorth.pdf" ) fn2 = file.path( odir, "plot_fishing_mortality_cfasouth.pdf" ) fn3 = file.path( odir, "plot_fishing_mortality_cfa4x.pdf" ) include_graphics(c(fn1, fn2, fn3) )

\@ref(fig:logisticFishingMortality)

## Stock status: Reference Points ... {.c}

```{r logistic-hcr, out.width='32%', echo=FALSE, fig.show='hold', fig.align='center', fig.cap = 'Reference Points (fishing mortality and modelled biomass) from Model 1, for N-ENS (left), S-ENS (middle), and 4X (right). The large yellow dot indicates most recent year and the 95\\% CI.' }
  odir = file.path( fishery_model_results, year.assessment, "logistic_discrete_historical" )
  fn1 = file.path( odir, 'plot_hcr_cfanorth.pdf' ) 
  fn2 = file.path( odir, 'plot_hcr_cfasouth.pdf' ) 
  fn3 = file.path( odir, 'plot_hcr_cfa4x.pdf' ) 
  include_graphics(c(fn1, fn2, fn3) )
#  \@ref(fig:logistic-hcr)

jae0/bio.snowcrab documentation built on Nov. 6, 2024, 10:10 p.m.