Random 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.
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.in
To 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++ (Debian 14.2.0-19) 14.2.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-3), silent = TRUE)
params.rndout <- gadgetutils::g3_optim(obj.rnd, params.rnd, trace = 0)
## initial value 406.457084
## iter 10 value 56.373102
## iter 20 value 10.773300
## iter 30 value 1.464858
## iter 40 value 0.884632
## iter 50 value 0.877144
## iter 60 value 0.867340
## iter 70 value 0.866218
## final value 0.865961
## converged
stopifnot(attr(params.rndout, "summary")$convergence)
Both models have ended up with similar solutions:
print(cbind(params.fixout$value, params.rndout$value))
## [,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 263.7658 619.9517
## fish_imm.init.2 2099.062 5306.031
## fish_imm.init.3 3815.14 6032.448
## fish_imm.init.4 113593.2 240268.1
## fish_imm.init.5 173688.4 45576.91
## fish_imm.M.1 0.9947778 0.7858061
## fish_imm.M.2 0.1001547 0.435409
## fish_imm.M.3 0.1007144 0.297616
## fish_imm.M.4 0.6172273 0.9999883
## fish_imm.M.5 0.3498054 0.284648
## init.F 0.8569412 0.6486456
## 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.127618e-09 1.175557e-06
## fish_imm.rec.proj.logar1.level 475.9299 881.6187
## fish_imm.rec.2000 484.58 903.6202
## fish_imm.rec.2001 482.8375 899.1764
## fish_imm.rec.2002 481.448 895.6371
## fish_imm.rec.2003 480.3393 892.8157
## fish_imm.rec.2004 479.4541 890.565
## fish_imm.rec.2005 478.7472 888.7685
## fish_imm.rec.2006 478.1824 887.3339
## fish_imm.rec.2007 477.7311 886.1879
## fish_imm.rec.2008 477.3703 885.2722
## fish_imm.rec.2009 477.0819 884.5403
## fish_imm.rec.2010 476.8512 883.9552
## fish_imm.rec.2011 476.6668 883.4874
## fish_imm.rec.2012 476.5194 883.1133
## fish_imm.rec.2013 476.4014 882.8142
## fish_imm.rec.2014 476.3071 882.575
## fish_imm.rec.2015 476.2316 882.3836
## fish_imm.rec.2016 476.1713 882.2306
## fish_imm.rec.2017 476.123 882.1082
## fish_imm.rec.2018 476.0844 882.0102
## fish_imm.rec.2019 476.0535 881.9319
## fish_imm.rec.2020 476.0287 881.8693
## fish_imm.rec.2021 476.009 881.8191
## fish_imm.rec.2022 475.9932 881.779
## fish_imm.rec.2023 475.9805 881.747
## fish_imm.rec.2024 475.9704 881.7213
## fish_imm.rec.2025 475.9623 881.7008
## fish_imm.rec.2026 475.9558 881.6844
## fish_imm.rec.2027 475.9506 881.6712
## fish_imm.rec.2028 475.9465 881.6607
## fish_imm.rec.2029 475.9432 881.6523
## fish_imm.rec.2030 475.9405 881.6456
## 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 0
Projections
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
## 17300000 36800000
Plot 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)
