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:
### Spawning & Random effects
library(gadget3)
library(dplyr)
actions <- list()
area_names <- g3_areas(c('IXa', 'IXb'))
# Create time definitions ####################
actions_time <- list(
g3a_time(
1980, 2000,
step_lengths = c(3L, 3L, 3L, 3L)),
NULL)
actions <- c(actions, actions_time)
# Create stock definition for fish ####################
st_imm <- g3_stock(c(species = "fish", 'imm'), seq(5L, 25L, 5)) |>
g3s_livesonareas(area_names["IXa"]) |>
g3s_age(1L, 5L)
st_mat <- g3_stock(c(species = "fish", 'mat'), seq(5L, 25L, 5)) |>
g3s_livesonareas(area_names["IXa"]) |>
g3s_age(3L, 10L)
stocks = list(imm = st_imm, mat = st_mat)
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)
# 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)
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)
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)
# Fleet data for f_surv #################################
# Landings data: For each year/step/area
expand.grid(year = 1980:2000, step = 2, area = 'IXa') |>
# Generate a random total landings by weight
mutate(weight = rnorm(n(), mean = 1000, sd = 100)) |>
# Assign result to landings_f_surv
identity() -> landings_f_surv
# Length distribution data: Generate 100 random samples in each year/step/area
expand.grid(year = 1980:2000, step = 2, area = 'IXa', length = rep(NA, 100)) |>
# Generate random lengths for these samples
mutate(length = rnorm(n(), mean = 50, sd = 20)) |>
# Save unagggregated data into ldist_f_surv.raw
identity() -> ldist_f_surv.raw
# Aggregate .raw data
ldist_f_surv.raw |>
# Group into length bins
group_by(
year = year,
step = step,
length = cut(length, breaks = c(seq(0, 80, 20), Inf), right = FALSE) ) |>
# Report count in each length bin
summarise(number = n(), .groups = 'keep') |>
# Save into ldist_f_surv
identity() -> ldist_f_surv
# Assume 5 * 5 samples in each year/step/area
expand.grid(year = 1980:2000, step = 2, area = 'IXa', age = rep(NA, 5), length = rep(NA, 5)) |>
# Generate random lengths/ages for these samples
mutate(length = rnorm(n(), mean = 50, sd = 20)) |>
# Generate random whole numbers for age
mutate(age = floor(runif(n(), min = 1, max = 5))) |>
# Group into length/age bins
group_by(
year = year,
step = step,
age = age,
length = cut(length, breaks = c(seq(0, 80, 20), Inf), right = FALSE) ) |>
# Report count in each length bin
summarise(number = n(), .groups = 'keep') ->
aldist_f_surv
# Create fleet definition for f_surv ####################
f_surv <- g3_fleet("f_surv") |> g3s_livesonareas(area_names["IXa"])
actions_f_surv <- list(
g3a_predate_fleet(
f_surv,
stocks,
suitabilities = g3_suitability_exponentiall50(by_stock = 'species'),
catchability_f = g3a_predate_catchability_totalfleet(
g3_timeareadata("landings_f_surv", landings_f_surv, "weight", areas = area_names))),
NULL)
actions_likelihood_f_surv <- list(
g3l_catchdistribution(
"ldist_f_surv",
obs_data = ldist_f_surv,
fleets = list(f_surv),
stocks = stocks,
function_f = g3l_distribution_sumofsquares(),
area_group = area_names,
report = TRUE,
nll_breakdown = TRUE),
g3l_catchdistribution(
"aldist_f_surv",
obs_data = aldist_f_surv,
fleets = list(f_surv),
stocks = stocks,
function_f = g3l_distribution_sumofsquares(),
area_group = area_names,
report = TRUE,
nll_breakdown = TRUE),
NULL)
actions <- c(actions, actions_f_surv, actions_likelihood_f_surv)
# Create abundance index for si_cpue ########################
# Generate random data
expand.grid(year = 1980:2000, step = 3, area = 'IXa') |>
# Fill in a weight column with total biomass for the year/step/area combination
mutate(weight = runif(n(), min = 10000, max = 100000)) ->
dist_si_cpue
actions_likelihood_si_cpue <- list(
g3l_abundancedistribution(
"dist_si_cpue",
dist_si_cpue,
stocks = stocks,
function_f = g3l_distribution_surveyindices_log(alpha = NULL, beta = 1),
area_group = area_names,
report = TRUE,
nll_breakdown = TRUE),
NULL)
actions <- c(actions, actions_likelihood_si_cpue)
# 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
# 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)