Introduction

SAMtool now supports MCMC sampling of the RCM using Stan called through the tmbstan R package. In this way, individual simulations within the operating model will be conditioned on individual draws of the MCMC chain.

Below is a schematic of the steps to setup the operating model using RCM and MCMC sampling. We will again use the Pacific cod example which sets up an age-structured model to mimic a delay-difference model, with knife-edge selectivity equal to knife-edge maturity.

MPD fit

The initial fit is generated by a function call to RCM() with a pre-made operating model (OM) and RCMdata object. This fit sets up the structure of the model, including any priors:

library(SAMtool)

#### Identify the maturity ogive to fix selectivity
data(pcod) 
mat_ogive <- pcod$OM@cpars$Mat_age[1, , 1]

#### Remove stochastic samples of M and steepness, will estimate from prior
OM <- pcod$OM
OM@cpars$M <- OM@cpars$h <- NULL
  
#### Update priors for M, h, and R0 (in addition to q)
prior <- pcod$prior
prior$R0 <- c(3, exp(1), exp(12))
prior$h <- c(0.7, 0.15)
prior$M <- c(0.5, 0.1)

#### MPD fit
MPD <- RCM(OM = OM, data = pcod$data, 
           condition = "catch", mean_fit = TRUE,
           selectivity = "free", s_selectivity = rep("SSB", ncol(pcod$data@Index)),
           start = list(vul_par = matrix(mat_ogive, length(mat_ogive), 1)),
           map = list(vul_par = matrix(NA, length(mat_ogive), 1),
                      log_early_rec_dev = rep(1, OM@maxage)),
           prior = prior)

The TMB model containing the MPD fit will be in MPD@mean_fit$obj.

MCMC

Once the MPD has been optimized, the posterior() function provides a convenient way to draw MCMC samples. This function is primarily a wrapper to tmbstan::tmbstan() which is also a wrapper to rstan::sampling(). Function arguments to posterior() configure the MCMC, e.g., the number of chains, iterations, thin rate, etc., and are passed on to these two functions.

MCMC <- posterior(MPD, chains = 2, iter = 2000, warmup = 1000, thin = 5)

The MCMC object is a stanfit object containing the parameters from the saved MCMC draws. Therefore, any methods and diagnostic functions that supports Stan models can be used with the MCMC object. These functions can be used to determine how whether the MCMC has converged.

#### Look at MCMC diagnostics (autocorrelation, wormplots, etc.)
sso <- shinystan::launch_shinystan(MCMC)

Update RCM

Now that we have a converged MCMC, we can update the operating model from a subset of these MCMC draws with RCMstan():

#### Identify the simulations and chains from which to grab parameters from the MCMC object
sims <- matrix(NA, OM@nsim, 2)
set.seed(234)
sims[, 1] <- sample(MCMC@sim$chains, OM@nsim, replace = TRUE)
sims[, 2] <- sample(MCMC@sim$n_save[1], OM@nsim)

#### Update RCM output and operating model
RCM_MCMC <- RCMstan(MPD, MCMC, sims)

RCMstan returns an updated RCModel and we can generate an updated markdown report:

plot(RCM_MCMC)

We can now run the MSE with the conditioned operating model.

MSE <- runMSE(RCM_MCMC@OM, ...)