Overview
This vignette demonstrates how to conduct Bayesian inference for the
CCSBT stock assessment model using the adnuts package and
the sparse inverse Hessian generated by RTMB. We illustrate
the setup, sampling procedure, and posterior diagnostics.
Setup and Optimization
# remotes::install_github("janoleko/RTMBdist")
# remotes::install_github("andrjohns/StanEstimators")
# remotes::install_github("Cole-Monnahan-NOAA/adnuts")
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Loading required package: RTMB
## Loading required package: RTMBdist
##
## Attaching package: 'RTMBdist'
##
## The following objects are masked from 'package:RTMB':
##
## dbeta, dnbinom2
##
## The following objects are masked from 'package:stats':
##
## dbeta, pnbinom, pt
##
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Loading required package: StanEstimators
theme_set(theme_bw())
attach(data_csv1)
lr <- data_labrep1
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,
removal_switch_f = c(0, 0, 0, 1, 0, 0), # 0=harvest rate, 1=direct removals
sel_min_age_f = c(2, 2, 2, 8, 6, 0, 2),
sel_max_age_f = c(17, 9, 17, 21, 25, 7, 17),
sel_end_f = c(1, 0, 1, 1, 1, 0, 1),
sel_LL1_yrs = c(1952, 1957, 1961, 1965, 1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997, 2001, 2006, 2007, 2008, 2011, 2014, 2017, 2020),
sel_LL2_yrs = c(1969, 2001, 2005, 2008, 2011, 2014, 2017, 2020),
sel_LL3_yrs = c(1954, 1961, 1965, 1969, 1970, 1971, 2005, 2006, 2007),
sel_LL4_yrs = c(1953),
sel_Ind_yrs = c(1976, 1995, 1997, 1999, 2002, 2004, 2006, 2008, 2010, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
sel_Aus_yrs = c(1952, 1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
sel_CPUE_yrs = c(1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997, 2001, 2006, 2007, 2008, 2011, 2014, 2017, 2020),
af_switch = 9,
lf_switch = 9, lf_minbin = c(1, 1, 1, 11),
cpue_switch = 1, cpue_a1 = 5, cpue_a2 = 17,
aerial_switch = 4, aerial_tau = 0.3,
troll_switch = 1,
pop_switch = 1,
hsp_switch = 1, hsp_false_negative = 0.7467647,
gt_switch = 1,
tag_switch = 1, tag_var_factor = 1.82
)
data <- get_data(data_in = data)
parameters <- get_parameters(data)
map <- get_map(data, parameters)
data$priors <- get_priors(parameters)
evaluate_priors(parameters, data$priors)## [1] 61.79457
obj <- MakeADFun(func = cmb(sbt_model, data), parameters = parameters, map = map)
bnd <- get_bounds(obj, parameters)
unique(names(obj$par))## [1] "par_log_B0" "par_log_m4" "par_log_m30"
## [4] "par_log_cpue_q" "par_log_troll_tau" "par_log_sel_1"
## [7] "par_log_sel_2" "par_log_sel_3" "par_log_sel_5"
## [10] "par_log_sel_6" "par_log_sel_7" "par_rdev_y"
obj$fn()## [1] 24545.05
control <- list(eval.max = 10000, iter.max = 10000)
opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, hessian = obj$he,
lower = bnd$lower, upper = bnd$upper, control = control)## outer mgc: 79.91711
## outer mgc: 287.5945
## outer mgc: 488.3391## Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
## NA/NaN function evaluation## outer mgc: 99.84999
## outer mgc: 6027.005## 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 evaluation## outer mgc: 10165.79
## outer mgc: 2230.533
## outer mgc: 150.7786
## outer mgc: 526.2794## Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
## NA/NaN function evaluation## outer mgc: 79.14928
## outer mgc: 612.4936## Warning in nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr, :
## NA/NaN function evaluation## outer mgc: 19.18882
## outer mgc: 5.25693
## outer mgc: 6.750307
## outer mgc: 10.69174
## outer mgc: 5.431459
## outer mgc: 5.431454
opt <- nlminb(start = opt$par, objective = obj$fn, gradient = obj$gr, hessian = obj$he,
lower = bnd$lower, upper = bnd$upper, control = control)## outer mgc: 5.431454
obj$par <- opt$parBayesian inferece
MCMC can be done using the sample_snuts function from
the adnuts package.
mcmc <- sample_snuts(
obj = obj, metric = "auto", init = "last.par.best",
# iter = 1000, warmup = 750, chains = 4, cores = 4,
num_samples = 10, num_warmup = 5, chains = 2, cores = 2,
control = list(adapt_delta = 0.99), globals = sbt_globals()
)## Optimizing...## Warning in .get_inputs(obj = obj, skip_optimization = skip_optimization, : Some
## standard errors estimated to be NaN, filling with dummy values so unit metric
## works. The 'mle' slot will be wrong so do not use it## Rebuilding RTMB obj without random effects...## Warning in sqrt(M): NaNs produced## unit metric selected b/c Qinv was not positive definite## log-posterior at inits=-23641.289; at conditional mode=-23641.289## Starting MCMC sampling...## Preparing parallel workspace...## Chain 2: Gradient evaluation took 0.010588 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 105.88 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: Iteration: 1 / 15 [ 6%] (Warmup)
## Chain 1: Gradient evaluation took 0.010734 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 107.34 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: Iteration: 1 / 15 [ 6%] (Warmup)
## Chain 1: Iteration: 6 / 15 [ 40%] (Sampling)
## Chain 1: Iteration: 15 / 15 [100%] (Sampling)
## Chain 1: Elapsed Time: 6.925 seconds (Warm-up)
## Chain 1: 0.186 seconds (Sampling)
## Chain 1: 7.111 seconds (Total)
## Chain 2: Iteration: 6 / 15 [ 40%] (Sampling)
## Chain 2: Iteration: 15 / 15 [100%] (Sampling)
## Chain 2: Elapsed Time: 7.048 seconds (Warm-up)
## Chain 2: 0.12 seconds (Sampling)
## Chain 2: 7.168 seconds (Total)## Registered S3 method overwritten by 'quantmod': method from as.zoo.data.frame zoo Attaching package: ‘RTMBdist’ The following objects are masked from ‘package:RTMB’: dbeta, dnbinom2 The following objects are masked from ‘package:stats’: dbeta, pnbinom, pt## 28 of 30 (93.33%) iterations ended with a divergence.
## These divergent transitions indicate that HMC is not fully able to explore the posterior distribution.
## Try increasing adapt_delta closer to 1.
## If this doesn't remove all divergences, try to reparameterize the model.##
##
## Model 'RTMB' has 1553 pars, and was fit using NUTS with a 'unit' metric
## 2 chain(s) of 15 total iterations (5 warmup) were used
## Average run time per chain was 7.14 seconds
## Minimum ESS=10 (50%), and maximum Rhat=1
## !! Warning: Signs of non-convergence found. Do not use for inference !!
## There were 20 divergences after warmup