Random effects & projections
Source:vignettes/articles/random-effect-recruitment.Rmd
random-effect-recruitment.RmdRandom effects in TMB run an inner optimisation problem for each run
of the model being optimised. To include random effects in a gadget3
model, you use the g3_param_project() series of
functions.
g3_param_project() works very much like
g3_parameterized(by_year = TRUE), introducing a time-vector
of parameters into a model. However, it also adds a likelihood component
describing a relationship between each element of the time-vector, which
is then used to pick values for projections.
For instance using
g3_param_project("M", g3_param_project_dnorm()) as the
definition of “M” in g3a_naturalmortality() will result in
the following parameters:
st <- g3_stock(c(species = "fish", "imm"), 1:10 * 10)
actions <- list(
g3a_time(2000, 2003, c(6,6)),
g3a_naturalmortality(st, g3a_naturalmortality_exp(
g3_param_project("M", g3_param_project_dnorm()) )))
model_cpp <- g3_to_tmb(actions)
attr(model_cpp, "parameter_template")[,c("switch", "random", "source")]## switch random source
## retro_years retro_years FALSE g3a_time
## proj_dnorm_M_weight proj_dnorm_M_weight FALSE g3_param_project
## M.proj.dnorm.mean M.proj.dnorm.mean FALSE g3_param_project_dnorm
## M.proj.dnorm.stddev M.proj.dnorm.stddev FALSE g3_param_project_dnorm
## M.2000.1 M.2000.1 TRUE g3_param_project
## M.2001.1 M.2001.1 TRUE g3_param_project
## M.2002.1 M.2002.1 TRUE g3_param_project
## M.2003.1 M.2003.1 TRUE g3_param_project
## M.2000.2 M.2000.2 TRUE g3_param_project
## M.2001.2 M.2001.2 TRUE g3_param_project
## M.2002.2 M.2002.2 TRUE g3_param_project
## M.2003.2 M.2003.2 TRUE g3_param_project
## project_years project_years FALSE g3a_time-
M.(year).(step): The time-vector of M values to use in historical timesteps -
M.proj.dnorm.mean/M.proj.dnorm.stddev: The parameters of the normal distribution we expectMto meet. These will…- Add the time-vector’s deviance from the ideal distribution to the model likelihood, so optimised solutions should fit the distribution.
- Be used to select future values of
Mwhen projecting
-
proj_dnorm_M_weight: The weighting of the likelihood component used above
Defining recruitment as a random effect
For a more complete example, we create a single-stock model, with a
g3l_abundancedistribution() likelihood component enforcing constant
abundance. We define the renewal rate using
g3_param_project() with
g3_param_project_logar1():
stocks <- list(
imm = g3_stock(c(species = "fish", "imm"), 1:10 * 10) |> g3s_age(1, 5) )
## Generate some abundnace data for ages 4 & 5, the stock should be stable, however the data is noisy
old.seed <- mget(c(".Random.seed"), envir = globalenv(), ifnotfound = list(NULL))[[1]]
set.seed(42)
dist_age4 <- expand.grid(year = 2000:2030, step = 1, age = 4, number = NA)
dist_age4$number <- rnorm(nrow(dist_age4), 1.75e7, 2e6)
dist_age5 <- expand.grid(year = 2000:2030, step = 1, age = 5, number = NA)
dist_age5$number <- rnorm(nrow(dist_age5), 3.75e7, 2e6)
set.seed(old.seed)
## Single stock model, using our generated data as an abundance index
actions <- list(
g3a_time(2000, 2030, step_lengths = c(6,6)),
g3a_age(stocks$imm),
g3a_initialconditions_normalcv(stocks$imm),
g3a_naturalmortality(stocks$imm),
g3a_growmature(stocks$imm, g3a_grow_impl_bbinom(maxlengthgroupgrowth = 4L)),
g3a_renewal_normalparam(
stocks$imm,
factor_f = g3_param_project(
"rec",
g3_param_project_logar1(),
scale = "scalar",
by_stock = stocks,
by_step = FALSE ),
run_step = 1 ),
g3l_abundancedistribution(
"dist_older",
rbind(dist_age4, dist_age5),
stocks = stocks,
# NB: Fix alpha/beta, abundance should match given values absolutely
function_f = g3l_distribution_surveyindices_log(alpha = 0, beta = 1),
report = TRUE,
nll_breakdown = TRUE),
NULL )
full_actions <- c(actions, list(
g3a_report_detail(actions),
g3a_report_history(actions, "__num$|__wgt$", out_prefix="dend_"), # NB: Late reporting
g3l_bounds_penalty(actions),
NULL ))
model_cpp <- g3_to_tmb(full_actions)Next, we define some initial parameters to use.
Note in particular that we are optimising the parameters for the distribution governing the time-vector. The aim here is roughly that the outer optimisation will decide on values here that meet the overall model, the inner optimistation will choose values that meet the distribution.
## Define baseline parameters / fixed-effects
attr(model_cpp, "parameter_template") |>
g3_init_val("*.init.scalar", 10, optimise = FALSE) |>
g3_init_val("*.init.#", 10, lower = 1e1, upper = 1e8) |>
g3_init_val("*.M.#", 0.5, lower = 0.1, upper = 1) |>
g3_init_val("init.F", 0.5, lower = 0.1, upper = 10) |>
g3_init_val("*_imm.Linf", 144.645) |>
g3_init_val("*.K", 0.3, lower = 0.04, upper = 1.2) |>
# NB: We don't optimise here since the model defines no length/weight structure, only counts
g3_init_val("*.t0", -0.8, optimise = FALSE) |>
g3_init_val("*.walpha", 0.01, optimise = FALSE) |>
g3_init_val("*.wbeta", 3, optimise = FALSE) |>
g3_init_val("*.rec.scalar", 10, optimise = FALSE) |>
g3_init_val("fish_imm.rec.proj.logar1.level", value = 2e3, lower = 1e1, upper = 5e5) |>
g3_init_val("fish_imm.rec.proj.logar1.stddev", value = 1, lower = 1e-15, upper = 1e4) |>
# NB: Relative weight is important, but g3_iterative isn't handling this for you yet
g3_init_val("proj_logar1_fish_imm_weight.proj_logar1_rec_weight", 1e0) |>
identity() -> params.inTo give something to compare against, optimise first with recruitment as a fixed-effect, which should still work, as the score for fixed & random effects is combined:
params.in |>
# NB: level = init value, to be internally consistent
g3_init_val("*.rec.#", value = 2e3, lower = 1e1, upper = 1e6, random = FALSE) |>
identity() -> params.fixed
obj.fix <- g3_tmb_adfun(model_cpp, params.fixed, silent = TRUE)## using C++ compiler: 'g++ (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0'
params.fixout <- gadgetutils::g3_optim(obj.fix, params.fixed, control = list(maxit = 2000, trace = 0))
stopifnot(attr(params.fixout, "summary")$convergence)Next, optimise the same model but using random effects for recruitment:
params.in |>
# NB: By setting random = TRUE, any previous bounds will be cleared
g3_init_val("*.rec.#", value = 2e3, random = TRUE) |>
identity() -> params.rnd
obj.rnd <- g3_tmb_adfun(model_cpp, params.rnd, inner.control = list(trace = 0, maxit = 1000, tol = 1e-7), silent = TRUE)
params.rndout <- gadgetutils::g3_optim(obj.rnd, params.rnd, trace = 0)## initial value 406.457084
## iter 10 value 56.336685
## iter 20 value 7.747918
## iter 30 value 1.200197
## iter 40 value 0.930482
## iter 50 value 0.894538
## iter 60 value 0.886271
## iter 70 value 0.881776
## final value 0.881775
## convergedBoth models have ended up with similar solutions:
## [,1] [,2]
## retro_years 0 0
## fish_imm.Linf 144.645 144.645
## fish_imm.K 0.3 0.3
## fish_imm.t0 -0.8 -0.8
## fish_imm.lencv 0.1 0.1
## fish_imm.init.scalar 10 10
## fish_imm.init.1 631.3684 1246.163
## fish_imm.init.2 9778.214 17214.09
## fish_imm.init.3 1135.981 9931.249
## fish_imm.init.4 8305.526 57666.91
## fish_imm.init.5 6539.932 95448.38
## fish_imm.M.1 0.9438641 0.9730494
## fish_imm.M.2 0.9980344 0.801599
## fish_imm.M.3 0.2635327 0.4103087
## fish_imm.M.4 0.6802546 0.6129333
## fish_imm.M.5 0.3469061 0.345459
## init.F 0.1791403 0.7280994
## recage 0 0
## fish_imm.walpha 0.01 0.01
## fish_imm.wbeta 3 3
## proj_logar1_fish_imm_weight.proj_logar1_rec_weight 1 1
## fish_imm.rec.proj.logar1.phi 0.8 0.8
## fish_imm.rec.proj.logar1.stddev 1.342989e-09 6.401645e-06
## fish_imm.rec.proj.logar1.level 1365.874 1289.954
## fish_imm.rec.2000 1397.306 1315.177
## fish_imm.rec.2001 1390.963 1310.094
## fish_imm.rec.2002 1385.908 1306.041
## fish_imm.rec.2003 1381.878 1302.807
## fish_imm.rec.2004 1378.662 1300.226
## fish_imm.rec.2005 1376.095 1298.165
## fish_imm.rec.2006 1374.045 1296.519
## fish_imm.rec.2007 1372.407 1295.203
## fish_imm.rec.2008 1371.098 1294.152
## fish_imm.rec.2009 1370.052 1293.311
## fish_imm.rec.2010 1369.215 1292.639
## fish_imm.rec.2011 1368.546 1292.101
## fish_imm.rec.2012 1368.012 1291.672
## fish_imm.rec.2013 1367.584 1291.328
## fish_imm.rec.2014 1367.242 1291.053
## fish_imm.rec.2015 1366.968 1290.833
## fish_imm.rec.2016 1366.749 1290.657
## fish_imm.rec.2017 1366.574 1290.516
## fish_imm.rec.2018 1366.434 1290.404
## fish_imm.rec.2019 1366.322 1290.314
## fish_imm.rec.2020 1366.233 1290.242
## fish_imm.rec.2021 1366.161 1290.184
## fish_imm.rec.2022 1366.104 1290.138
## fish_imm.rec.2023 1366.058 1290.101
## fish_imm.rec.2024 1366.021 1290.072
## fish_imm.rec.2025 1365.992 1290.048
## fish_imm.rec.2026 1365.968 1290.029
## fish_imm.rec.2027 1365.95 1290.014
## fish_imm.rec.2028 1365.935 1290.002
## fish_imm.rec.2029 1365.923 1289.992
## fish_imm.rec.2030 1365.913 1289.985
## report_detail 1 1
## fish_imm.bbin 0 0
## fish_imm.rec.scalar 10 10
## fish_imm.rec.sd 10 10
## adist_surveyindices_log_dist_older_weight 1 1
## project_years 0 0Projections
Given a parameterised model, projections can be performed by varying
the project_years parameter. Whilst you can use the
objective function above for this, it’s simpler to do so with a function
designed for single runs, produced either by g3_tmb_fn or
g3_tmb_r. We will use the latter here:
fn <- g3_tmb_fn(model_cpp)The resulting fn function then accepts parameter tables,
which we can modify to project into the future. Our model settles at the
point we expect:
abund <- g3_array_agg(
fn(params.rndout |> g3_init_val("project_years", 40))$dend_fish_imm__num,
c('age', 'year'),
age = c(4, 5),
year = 2040:2060,
step = 1)
print(signif(rowMeans(abund), 3))## age4 age5
## 17400000 36800000Plot abundance for years 4 & 5, with 40 years projection, using level/stddev chosen by optimiser:
r <- fn(params.rndout |> g3_init_val("project_years", 40))
g3_array_plot(t(g3_array_agg(r$dend_fish_imm__num, c('age', 'year'))))
Successive calls will choose different values for the projections. Plot 20 iterations of projecting 40 years:
proj_iters <- do.call(cbind, lapply(1:20, function (i) {
r <- fn(params.rndout |> g3_init_val("project_years", 40))
# Generate report for age 4 abundance
proj <- t(g3_array_agg(
r$dend_fish_imm__num,
c('age', 'year'),
year = 2025:2060, # NB: Skip the first years when the model is settling
age = 4,
step = 1 ))
dimnames(proj)[[2]] <- paste0("proj ", i)
return(proj)
}))
g3_array_plot(proj_iters)