Load inputs
Below the required R libraries are loaded. The sbt TMB
model is dynamically loaded along with several R functions using
library(sbt). The kableExtra package is used
for generating tables of inputs. The code
theme_set(theme_bw()) alters the plot aesthetics.
Below are example inputs that are loaded from the sbt
package. Normally the user would import these inputs (e.g., from
.csv files using the read_csv function).
names(data_par1)
#> [1] "ln_B0" "lnq" "lnqhsp"
#> [4] "deltalnq86" "steep" "sigma_r"
#> [7] "sigma_cpue" "m4" "m0"
#> [10] "m10" "m30" "cpue_omega"
#> [13] "psi" "Reps" "ln_bdiff"
#> [16] "ln_sel_aerial" "tau_aerial" "tau_troll"
#> [19] "tag_H_factor" "par_sels_init_i" "par_sels_change_i"
#> [22] "par_log_hstar_i"
names(data_labrep1)
#> [1] "scenario_number" "n" "years"
#> [4] "ObjF" "lnlike" "penal"
#> [7] "Hhigh" "H2003" "H2004"
#> [10] "sigma.cpue" "tag.var.factor" "gtOD.factor"
#> [13] "tau.aerial" "tau.troll" "sigma.r"
#> [16] "rec_AC" "AC_penalty" "q_AC"
#> [19] "res.stats" "B0" "R0"
#> [22] "alpha" "beta" "rho"
#> [25] "steep" "depletion" "psi"
#> [28] "lnq" "qhsp" "qgt"
#> [31] "M" "omega" "TOTbio"
#> [34] "B10_plus" "Sbio" "phi"
#> [37] "phsp" "Recruitment" "Rdev"
#> [40] "Reps" "ages" "yrs.catch"
#> [43] "Ns" "Fs" "sel"
#> [46] "yrs.cpue" "cpue" "cpue.pred"
#> [49] "cpue.resid" "Aerial.sw" "Aerial.Surv"
#> [52] "Troll.sw" "Troll.Index" "ages.1"
#> [55] "age.obs.1" "age.pred_1" "age.res_1"
#> [58] "ages.2" "age.obs.2" "age.pred_2"
#> [61] "age.res_2" "lengths" "len.obs"
#> [64] "len.pred" "len.res" "catch.pred"
#> [67] "catch.weight.pred" "specs" "Nsamp"
#> [70] "sel.change.sd" "sel.smooth.sd" "tag.obs"
#> [73] "tag.pred" "tag.res" "pooled.tag.obs"
#> [76] "pooled.tag.pred" "pooled.1yeartag.obs" "pooled.1yeartag.pred"
#> [79] "spweights.age" "weights.ageLL1" "weights.ageLL2"
#> [82] "weights.ageSurf" "Bage" "Fage"
#> [85] "MSY" "Bmsy" "TBmsy"
#> [88] "Fmsy_a215" "F_a215" "Fmsy_a"
#> [91] "Bmsy_a" "sB0" "sB0.sd"
#> [94] "sB1980" "sB1980.sd" "sBfinal"
#> [97] "sBfinal.sd" "dep1" "dep1.sd"
#> [100] "dep2"Define the data list. The get_data function sets up some
additional inputs that are required before the data list can be passed
to MakeADFun.
Data <- list(last_yr = 2022, age_increase_M = 25,
length_m50 = 150, length_m95 = 180,
catch_UR_on = 0, catch_surf_case = 1, catch_LL1_case = 1,
scenarios_surf = scenarios_surface, scenarios_LL1 = scenarios_LL1,
sel_min_age_f = c(2, 2, 2, 8, 6, 0),
sel_max_age_f = c(17, 9, 17, 22, 25, 7),
sel_end_f = c(1, 0, 1, 1, 1, 0),
sel_change_sd_fy = t(as.matrix(sel_change_sd[,-1])),
sel_smooth_sd_f = data_labrep1$sel.smooth.sd,
removal_switch = c(0, 0, 0, 0, 0, 0), # 0=standard removals, 1=direct removals
pop_switch = 1,
hsp_switch = 1, hsp_false_negative = 0.7467647,
gt_switch = 1,
cpue_switch = 1, cpue_a1 = 5, cpue_a2 = 17,
aerial_switch = 4, aerial_tau = 0.3,
troll_switch = 0,
af_switch = 3,
lf_switch = 3, lf_minbin = c(1, 1, 1, 11),
tag_switch = 1, tag_var_factor = 1.82
)
Data <- get_data(data_in = Data)Input tables
Below are examples of some of the model inputs in tabular form.
Mean length (cm) by year, season, and age. Only a subset of years and ages 1-10 are shown:
kbl(bind_rows(head(length_mean), tail(length_mean))[,1:13]) %>%
add_header_above(c(" " = 2, "Age" = 11))| Year | Season | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1931 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 1932 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 1933 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 1934 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 1935 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 1936 | 1 | 45 | 57.1 | 74.1 | 88.8 | 101.8 | 113.9 | 124.3 | 133.3 | 141.0 | 147.6 | 153.3 |
| 2017 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
| 2018 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
| 2019 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
| 2020 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
| 2021 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
| 2022 | 2 | 50 | 68.9 | 91.6 | 105.5 | 117.3 | 127.2 | 135.6 | 142.7 | 148.7 | 153.8 | 158.1 |
The parent offspring pair (POP) data set:
| Cohort | CaptureYear | CaptureAge | NPOPs | Comps |
|---|---|---|---|---|
| 2004 | 2021 | 30 | 0 | 6625 |
| 2005 | 2021 | 30 | 0 | 6780 |
| 2006 | 2021 | 30 | 0 | 6735 |
| 2007 | 2021 | 30 | 0 | 6575 |
| 2008 | 2021 | 30 | 0 | 4815 |
| 2009 | 2021 | 30 | 0 | 4380 |
| 2010 | 2021 | 30 | 0 | 4515 |
| 2011 | 2021 | 30 | 0 | 4495 |
| 2012 | 2021 | 30 | 1 | 4765 |
| 2013 | 2021 | 30 | 0 | 4270 |
| 2014 | 2021 | 30 | 0 | 4740 |
| 2015 | 2021 | 30 | 0 | 3780 |
| 2016 | 2021 | 30 | 0 | 7245 |
| 2017 | 2021 | 30 | 0 | 7560 |
| 2018 | 2021 | 30 | 0 | 6920 |
The unaccounted catch (tonnes) data set:
| Year | LL1 | LL2 | LL3 | LL4 | Indonesia | Australia |
|---|---|---|---|---|---|---|
| 2008 | 72 | 0 | 0 | 0 | 0 | 0 |
| 2009 | 152 | 0 | 0 | 0 | 0 | 0 |
| 2010 | 271 | 0 | 0 | 0 | 0 | 0 |
| 2011 | 151 | 0 | 0 | 0 | 0 | 0 |
| 2012 | 275 | 0 | 0 | 0 | 0 | 0 |
| 2013 | 432 | 0 | 0 | 0 | 0 | 0 |
| 2014 | 121 | 0 | 0 | 0 | 0 | 0 |
| 2015 | 326 | 0 | 0 | 0 | 0 | 0 |
| 2016 | 756 | 0 | 0 | 0 | 0 | 0 |
| 2017 | 984 | 0 | 0 | 0 | 0 | 0 |
| 2018 | 1511 | 0 | 0 | 0 | 0 | 0 |
| 2019 | 1155 | 0 | 0 | 0 | 0 | 0 |
| 2020 | 1160 | 0 | 0 | 0 | 0 | 0 |
| 2021 | 1160 | 0 | 0 | 0 | 0 | 0 |
| 2022 | 1160 | 0 | 0 | 0 | 0 | 0 |
Input plots
plot_length_at_age(data = Data)
Length (cm) at age during each season.
plot_weight_at_age(data = Data)
Weight (kg) at age for each fishery.
Model setup
Parameters
Define the parameter list:
Params <- list(par_log_B0 = data_par1$ln_B0,
par_log_psi = log(data_par1$psi),
par_log_m0 = log(data_par1$m0),
par_log_m4 = log(data_par1$m4),
par_log_m10 = log(data_par1$m10),
par_log_m30 = log(data_par1$m30),
par_log_h = log(data_par1$steep),
par_log_sigma_r = log(data_labrep1$sigma.r),
par_log_cpue_q = data_par1$lnq,
par_log_cpue_sigma = log(data_par1$sigma_cpue),
par_log_cpue_omega = log(data_par1$cpue_omega),
par_log_aerial_tau = log(data_par1$tau_aerial),
par_log_aerial_sel = data_par1$ln_sel_aerial,
par_log_troll_tau = log(data_par1$tau_troll),
par_log_hsp_q = data_par1$lnqhsp,
par_log_tag_H_factor = log(data_par1$tag_H_factor),
par_log_af_alpha = c(1, 1),
par_log_lf_alpha = c(1, 1, 1, 1),
par_rdev_y = data_par1$Reps,
par_sels_init_i = data_par1$par_sels_init_i,
par_sels_change_i = data_par1$par_sels_change_i
)Use TMB’s Map option to turn parameters on/off:
Map <- list()
Map[["par_log_psi"]] <- factor(NA)
Map[["par_log_m0"]] <- factor(NA)
Map[["par_log_m10"]] <- factor(NA)
Map[["par_log_h"]] <- factor(NA)
Map[["par_log_sigma_r"]] <- factor(NA)
Map[["par_log_cpue_sigma"]] <- factor(NA)
Map[["par_log_cpue_omega"]] <- factor(NA)
Map[["par_log_aerial_tau"]] <- factor(NA)
Map[["par_log_aerial_sel"]] <- factor(rep(NA, 2))
Map[["par_log_troll_tau"]] <- factor(NA)
Map[["par_log_hsp_q"]] <- factor(NA)
Map[["par_log_tag_H_factor"]] <- factor(NA)
Map[["par_log_af_alpha"]] <- factor(rep(NA, 2))
Map[["par_log_lf_alpha"]] <- factor(rep(NA, 4))Specify the random effects:
Random <- c()Normally, recruitment deviates (i.e., par_rdev_y) should
be treated as random effects, but they are left as fixed effects for
comparability.
Create the AD object using TMBs MakeADFun function:
Optimisation
List of parameters that are on:
unique(names(obj$par))
#> [1] "par_log_B0" "par_log_m4" "par_log_m30"
#> [4] "par_log_cpue_q" "par_rdev_y" "par_sels_init_i"
#> [7] "par_sels_change_i"The objective function value:
obj$fn(obj$par)
#> [1] 27328.36Load the default parameter bounds:
bnd <- get_bounds(obj = obj)Optimise using the nlminb function:
opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr,
lower = bnd$lower, upper = bnd$upper)
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluation
#> Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
#> NA/NaN function evaluationCalculate standard deviations of all model parameters, including non
linear functions of random effects and parameters specified through the
ADREPORT() macro from the user template:
Check that the Hessian is positive definite:
print(Report$pdHess)
#> [1] FALSETry the analytical Hessian instead:
# he <- obj$he()
# he_inv <- solve(he)
# he_ch <- chol(he)
# ev <- eigen(he)
# range(ev$values)Plot outputs
plot_natural_mortality(data = Data, object = obj)
plot_initial_numbers(data = Data, object = obj)
plot_catch(data = Data, object = obj)
#> [1] "The maximum catch difference was: 1.78220943780616e-07"
plot_hrate(data = Data, object = obj, years = 1990:2010)
# yrs <- c(1952, 1953, 1954, seq(1955, 2020, 5), 2021, 2022)
yrs <- c(1952:2022)
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "LL1")
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "LL2")
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "LL3")
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "LL4")
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "Indonesian")
plot_selectivity(data = Data, object = obj, years = yrs, fisheries = "Australian")
plot_recruitment(data = Data, object = obj)
plot_biomass_spawning(data = Data, object = obj)
