```{css, echo=FALSE} / https://bookdown.org/yihui/rmarkdown-cookbook/html-css.html / .modelScript pre, .modelScript pre.sourceCode code { background-color: #f0f8ff !important; }
```r library(unittest) # Redirect ok() output to stderr options(unittest.output = stderr()) if (nzchar(Sys.getenv('G3_TEST_TMB'))) options(gadget3.tmb.work_dir = gadget3:::vignette_base_dir('work_dir')) library(gadget3) set.seed(123)
This vignette walks through a script that will generate a gadget3 model, explaining concepts along the way.
Code blocks that make up the gadget3 script are marked in blue, like this:
When combined they will form a full model, see the appendix for the entire script.
As mentioned before in vignette('introduction-single-stock')
,
gadget3 stock objects do not have to correspond 1:1 with a species.
We can have multiple stock objects representing the same species in a different stage in their life-cycle,
most commonly mature
and immature
versions, male
and female
versions, or all 4.
The set up is much the same as before, but major differences will be highlighted.
Initial setup & time-keeping is identical:
We define 2 stocks instead of one, and a list()
containing both for convenience:
Notice that:
"fish_imm"
.Stock actions now need to include interactions between immature & mature:
actions_st_imm
and actions_st_imm
are largely similar to our actions_fish
from the previous model, but:
maturity_f
to g3a_growmature()
to move individuals to the mature stock.
The rate of maturity is coupled to growth, which is why g3a_growmature()
does both at the same time.g3a_age()
can also move individuals to the mature stock.
This will happen if an immature fish ages beyond the final age bin (5 in our case).
At that point it matures "by default".g3a_renewal_normalcv()
, as there is no recruitment directly into the mature stock.There is very little difference defining a fleet for a multiple stock model vs. single stocks.
To define a fleet, we need to introduce historical data into the model. In our case we will generate random data to use later:
For more information on how this works, see vignette("incorporating-observation-data")
.
As well as the data we defined last time, we also define a maturity stage dataset. We represent the maturity stage data by adding a "stock" column, when we then compare to model gadget3 will breakdown predicted catches by stock.
matp_f_surv[sample(nrow(matp_f_surv), 10),]
Also note we don't provide the full name of a stock, this both makes code re-use easier and means that we would combine multiple stocks if there is for example a "fish_imm_m" and "fish_imm_f".
Our fleet is defined with the same set of actions as the single-species model:
There are 2 differences to before:
stocks
, not an individual stock.
Without additional changes, the 2 are treated as a single combined stock.g3_suitability_exponentiall50(by_stock = 'species')
,
instructing it to have a single alpha
& l50
parameter for both our stocks,
as they have the same species name.g3l_catchdistribution()
as before,
the difference is in the columns of our input data.The by_stock
parameter is passed through to g3_parameterized()
.
We can see the result of setting this in the parameter template.
If by_stock = TRUE
(the default) then we get parameters for both fish_imm.f_surv.l50
& fish_mat.f_surv.l50
:
suppressWarnings({ # NB: The model is incomplete, fish_imm__num isn't defined simple_model <- g3_to_r(list(g3a_time(1990, 2023), g3a_predate_fleet( f_surv, stocks, suitabilities = g3_suitability_exponentiall50(by_stock = TRUE), catchability_f = g3a_predate_catchability_totalfleet(1) ))) }) names(attr(simple_model, "parameter_template"))
ok(ut_cmp_identical(sort(names(attr(simple_model, "parameter_template")), method = "radix"), c( "fish_imm.f_surv.alpha", "fish_imm.f_surv.l50", "fish_mat.f_surv.alpha", "fish_mat.f_surv.l50", "project_years", "retro_years", NULL)), "Params for simple_model, by_stock = TRUE")
If by_stock = 'species'
, then there is a single, shared fish.f_surv.l50
parameter:
suppressWarnings({ # NB: The model is incomplete, fish_imm__num isn't defined simple_model <- g3_to_r(list(g3a_time(1990, 2023), g3a_predate_fleet( f_surv, stocks, suitabilities = g3_suitability_exponentiall50(by_stock = 'species'), catchability_f = g3a_predate_catchability_totalfleet(1) ))) names(attr(simple_model, "parameter_template")) })
ok(ut_cmp_identical(sort(names(attr(simple_model, "parameter_template")), method = "radix"), c( "fish.f_surv.alpha", "fish.f_surv.l50", "project_years", "retro_years", NULL)), "Params for simple_model, by_stock = 'species'")
The by_stock
parameter is just a convenient shortcut to change the default settings,
if we specify g3_parameterized()
ourselves we can change the parameterization in other ways,
for example by_year = TRUE
gives us per-year l50
:
suppressWarnings({ # NB: The model is incomplete, fish_imm__num isn't defined simple_model <- g3_to_r(list(g3a_time(1990, 2023), g3a_predate_fleet( f_surv, stocks, suitabilities = g3_suitability_exponentiall50( l50 = g3_parameterized("l50", by_stock = 'species', by_predator = TRUE, by_year = TRUE)), catchability_f = g3a_predate_catchability_totalfleet(1) ))) names(attr(simple_model, "parameter_template")) })
ok(ut_cmp_identical(sort(names(attr(simple_model, "parameter_template")), method = "radix"), c( paste0("fish.f_surv.l50.", 1990:2023), "fish_imm.f_surv.alpha", "fish_mat.f_surv.alpha", "project_years", "retro_years", NULL)), "Params for simple_model, by_year = TRUE")
See ?vignette('model-customisation')
for more.
If we had data on the distribution of mature vs. immature, our observation data could contain a stock
column with fish_imm
or fish_mat
. See vignette("incorporating-observation-data")
.
Again, further fleets can be added by repeating the code above.
Survey indices should be handed the full stocks
list, instead of a single stock,
but are otherwise the same as before:
At this point, we are ready to convert our model into code:
Now we should be configuring parameters based on the template.
Thanks to using wildcards in the g3_init_val()
calls,
a lot of the parameter settings will work regardless of the model being single- or multi-stock,
so we don't need to change the initial values from the previous model:
Finally we are ready for optimisation runs.
g3_tmb_adfun()
is a wrapper around TMB::MakeADFun()
and TMB::compile
, producing a TMB objective function.
gadgetutils::g3_iterative()
then optimises based on iterative reweighting
Once this has finished, we can view the output using gadgetplots::gadget_plots()
.
Once finished, you can view the output in your web browser:
For convenience, here is all the sections of the model script above joined together:
```{js, echo=FALSE} document.write( '
' + Array.from(document.querySelectorAll('.modelScript pre')).map((x) => x.innerHTML).join("\n\n") + '
```r gadget3:::vignette_test_output("multiple-substocks", model_code, params.out)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.