NB: This vignette is a work in progress
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:
### Spawning & Random effects
We only show the blocks that are relevant to implementing random effects, see the appendix for the entire script.
For an explanation of other parts of the script, see previous vignettes.
In vignette("introduction-single-stock")
we introduced g3a_renewal_normalparam()
as a means of adding recruitment into the model.
An alternative action is g3a_spawn()
. The main difference here is the existing state of the stock is used to inform the rate of recruitment.
To improve our model fit, we shall define the rate of spawning using a random walk.
Firstly, we define our immature stock. There’s no real difference here, apart from g3a_renewal_normalparam()
no longer being present:
actions_st_imm <- list(
g3a_growmature(st_imm,
g3a_grow_impl_bbinom(
maxlengthgroupgrowth = 4L ),
# Add maturation
maturity_f = g3a_mature_continuous(),
output_stocks = list(st_mat),
transition_f = ~TRUE ),
g3a_naturalmortality(st_imm),
g3a_initialconditions_normalcv(st_imm),
# NB: g3a_renewal_normalparam() no longer here
g3a_age(st_imm, output_stocks = list(st_mat)),
NULL)
actions <- c(actions, actions_st_imm)
Next, we need to define the set of parameters that will constitute our random walk:
# Configure our random walk for mu
st_spawn_mu <- g3_parameterized('spawn_mu',
by_year = TRUE,
scale = g3_parameterized('spawn_mu_seasonal', by_step = TRUE, value = 1, optimise = FALSE),
random = TRUE)
As we saw in vignette('multiple-substocks')
, by_year
will give us a series of parameters for spawn_mu
, one per year. random = TRUE
will ensure that these parameters have TRUE in the random column in the parameter template.
The scale
parameter adds seasonal variation parameters, spawn_mu_seasonal.1
..spawn_mu_seasonal.4
, which can be used to control when spawning happens.
We store this as a separate variable as we need to use it in 2 places. First, when we configure our spawning action:
actions_st_mat <- list(
g3a_growmature(st_mat,
g3a_grow_impl_bbinom(
maxlengthgroupgrowth = 4L )),
g3a_naturalmortality(st_mat),
g3a_initialconditions_normalcv(st_mat),
g3a_spawn(
st_mat,
recruitment_f = g3a_spawn_recruitment_bevertonholt(
mu = st_spawn_mu,
lambda = g3_parameterized("spawn_lambda", by_stock = TRUE) ),
proportion_f = g3_suitability_exponentiall50(),
output_stocks = list(st_imm) ),
g3a_age(st_mat),
NULL)
actions <- c(actions, actions_st_mat)
Mature stocks spawn, their output going into the immature stock. We use suitability functions to define proportion_f
, the proportion of stock ready to spawn, by_step
meaning this proportion is seasonal. Our random walk is used as the rate of recruitment, mu.
Next we use a likelihood action to constrain the values in the random walk during optimisation:
actions_likelihood_st <- list(
g3l_understocking(stocks, nll_breakdown = TRUE),
g3l_random_walk('rwalk_st_spawn_mu',
param_f = st_spawn_mu,
sigma_f = g3_parameterized('spawn_mu.sigma', value = 5, optimise = FALSE),
log_f = FALSE ),
g3l_random_dnorm('rwalk_st_spawn_mu',
param_f = st_spawn_mu,
mean_f = g3_parameterized('spawn_mu.mean', value = 78, optimise = FALSE),
sigma_f = g3_parameterized('spawn_mu.sigma', value = 1, optimise = FALSE),
log_f = FALSE ),
NULL)
actions <- c(actions, actions_likelihood_st)
For random effects to be useful, they have to be constrained to sensible bounds. There are 2 likelihood actions that help with this, g3l_random_dnorm()
& g3l_random_walk()
. The former penalises random values that stray outside given normal distribution, the latter constrains along a random walk; the current year’s value should not vary significantly from the previous.
At this point, we are ready to convert our model into code:
# Create model objective function ####################
model_code <- g3_to_tmb(c(actions, list(
g3a_report_detail(actions),
g3l_bounds_penalty(actions) )))
# Guess l50 / linf based on stock sizes
estimate_l50 <- g3_stock_def(st_imm, "midlen")[[length(g3_stock_def(st_imm, "midlen")) / 2]]
estimate_linf <- max(g3_stock_def(st_imm, "midlen"))
estimate_t0 <- g3_stock_def(st_imm, "minage") - 0.8
attr(model_code, "parameter_template") |>
g3_init_val("*.rec|init.scalar", 10, lower = 0.001, upper = 200) |>
g3_init_val("*.init.#", 10, lower = 0.001, upper = 200) |>
g3_init_val("*.M.#", 0.15, lower = 0.001, upper = 1) |>
g3_init_val("init.F", 0.5, lower = 0.1, upper = 1) |>
g3_init_val("*.Linf", estimate_linf, spread = 0.2) |>
g3_init_val("*.K", 0.3, lower = 0.04, upper = 1.2) |>
g3_init_val("*.t0", estimate_t0, spread = 2) |>
g3_init_val("*.walpha", 0.01, optimise = FALSE) |>
g3_init_val("*.wbeta", 3, optimise = FALSE) |>
g3_init_val("*.*.alpha", 0.07, lower = 0.01, upper = 0.2) |>
g3_init_val("*.*.l50", estimate_l50, spread = 0.25) |>
g3_init_val("*.bbin", 100, lower = 1e-05, upper = 1000) |>
g3_init_val("spawn_mu.#", value = 78) |>
#g3_init_val("spawn_mu.1980", value = 78, random = FALSE) |>
#g3_init_val("spawn_lambda", value = 1e-6, optimise = TRUE) |>
identity() -> params.in
# Uncomment this to temporarily disable random effects
#params.in[params.in$random, 'optimise'] <- TRUE
#params.in[params.in$random, 'random'] <- FALSE
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
# Optimise model ################################
obj.fn <- g3_tmb_adfun(model_code, params.in, inner.control = list(trace = 3, maxit = 100))
#obj.fn$env$tracepar <- TRUE
#obj.fn$env$tracemgc <- TRUE
out <- optim(par = obj.fn$par, fn = obj.fn$fn, gr = obj.fn$gr, method = 'BFGS', control = list(
maxit = 1000,
trace = 1,
reltol = .Machine$double.eps^2 ))
params.out <- g3_tmb_relist(params.in, out$par)
fit <- gadgetutils::g3_fit(model_code, params.out)
For convenience, here is all the sections of the model script above joined together: