Skip to content

3: Ricker Model Forms in BUGS & JAGS

Jump to bottom
SOLV-Code edited this page Jan 10, 2025 · 1 revision

Introduction

The RapidRicker package currently includes Bayesian implementations for 3 forms of the Ricker model:

  • Ricker (Basic Ricker Model): density-dependent productivity, but relationship is assumed to be stable over time (i.e. R/S changes with spawner abundance, but for a given spawner abundance the expected recruits are the same in 1980 and in 2005. -> get 1 alpha parameter (productivity) and 1 beta parameter (capacity)
  • RickerAR1 (Ricker Model with Lag1 autoregression): checks for an underlying pattern in the productivity residuals, then removes the pattern based on the assumption that there is a stable underlying relationship, before estimating the parameters -> get 1 alpha parameter (productivity) and 1 beta parameter (capacity)
  • RickerKalman (Ricker model with time-varying productivity): checks for an underlying pattern in the productivity residuals, then estimates a year-specific productivity parameter -> get 1 alpha parameter for each brood year (productivity) and 1 beta parameter that is fixed over time (capacity)

BUGS/JAGS Code

Concepts and Definitions

  • BUGS and JAGS are two different programs for doing the random sampling.
  • Both use very similar (and mostly interchangeable) model descriptions (i.e. if your model runs in BUGS, it should run in JAGS, and the other way around).
  • BUGS/JAGS models mostly look like R code, but there are some key differences in the functions, especially when specifying probability distributions ( LeBauer, Dietze, and Bolker 2013). In R, you sample 50 values from a normal distribution with mean = 0 and sd = 2 as rnorm(50,0,2), because the R function takes the arguments dnorm(x,mean,sd). In BUGS/JAGS, the same sample is generated by ~dnorm(0,0.25) because the BUGS/JAGS function takes arguments dnorm(mean,precision), where the precision parameter tau = (1/sd)^2.

Basic Ricker

The Basic Ricker JAGS code (available here) was originally developed by Catherine Michielsens (PSC) and has been widely used in Pacific salmon analyses (e.g. Fraser River sockeye forecast, Fraser River sockeye Recovery Potential Assessment, various transboundary escapement goal analyses, INCL LINKS) with various case-specific modifications. The current version included in the RapidRicker package uses the notation from previous northern transboundary JAGS implementations, combined with a recent addition by Ann-Marie Huang (DFO) that implemented an upper constraint on the capacity prior (see below).

The key feature of the JAGS code is the basic Ricker equation, which in JAGS has a few distinct steps:

  • Observed Recruits R_Obs are a function of predicted recruits plus noise
  • Predicted log recruits logR are a function of predicted recruits/spawner (in log space) times the log of the brood year spawner abundance
  • Predicted recruits/spawner are a function of the productivity parameter ln.alpha and the capacity paramter beta.

In JAGS code that looks like:

      R_Obs[i] ~ dlnorm(logR[i],tau_R)          
      logR[i] <- RS[i] +log(S[i])               
      RS[i] <- ln.alpha - beta * S[i]

Ricker AR1

The Ricker AR1 JAGS code (available here) has been expanded for AR1 from the Basic Ricker code described above based on Eq21 and 22 of Fleischman and Evenson (2010; ADFG FMS10-04). An example of previous use is this Taku Coho Benchmarks paper.

The key feature is an extension of the basic Ricker equation with a year-specific residual that depends on the residual from 1 year earlier and a scaling parameter phi. This way "good" years (positive residuals) will tend to follow "good" years, and "bad" years tend to follow "bad" years, capturing an underlying pattern over time, but excluding it from the estimates of alpha and beta.

RS[i] <- ln.alpha - beta * S[i] + phi * log.resid[i-1]

Ricker Kalman

The Ricker Kalman JAGS code (available here) has been adapted by Ann-Marie Huang (DFO) from the Basic Ricker code described above based on SOURCE. An example of previous use is the Recovery Potential Assessment of 9 Fraser Sockeye Conservation Units.

The key feature is an extension of the basic Ricker equation with a year-specific productivity parameter alpha that depends on the productivity from 1 year earlier and an incremental change parameter w. The precision parameter tauw specifies how big the change in productivity can be from 1 year to the next.

 RS[i] <- ln.alpha[i] - beta * S[i]   

where

ln.alpha[i] <- ln.alpha[i-1] + w[i]
w[i]~ dnorm(0,tauw)

Arguments in the calcMCMCRickerBM() Function

Basic Inputs and Settings

Variable Definition Sample Value
sr_obj spawner-recruit data data frame with (at least) the variables Year, Spn and Rec. Others are ignored.
sr.scale scale paramter for the SR data an integer value used to rescale the Spn and Rec variables in sr_obj, prior to the MCMC fit, default = 10^6 (i.e. convert to millions). NOTE: If sr.scale is different from 1, then the benchmark estimates are scaled back up, but the MCMC estimates of alpha and beta will be in different units then the alpha and beta estimates from the deterministic fit.
model.type type of model being fitted "RickerBasic" of "RickerAR1" or "RickerKalman"

To be incorporated:

  • model.file = "BUILT_IN_MODEL_Ricker_BUGS.txt",

  • min.obs=15, ,

      			output = "short",

  			out.path = "MCMC_Out",

  			out.label = "MCMC",

  			mcmc.seed = "default",

  			tracing = FALSE

MCMC Settings

Variable Definition Sample Value
n.chains number of separate sequences of MCMC samples 2
n.burnin number of initial samples to discard 20,000
n.thin number of subsamples 60
n.samples TOTAL number of samples per chain 50,000

These are included as a list in the function call:

mcmc.settings = list(n.chains=2, n.burnin=20000, n.thin=60,n.samples=50000)

Using the sample settings listed above:

  • For each chain:
    • 50k total samples - 20k burnin = 30k used samples
    • 30k used samples / 60 thinning = 500 retained samples
  • 500 samples / chain * 2 chains = 1,000 samples in the output
  • to increase the burn.in, you need to increase BOTH the burnin and the total sample size. To increase the sample size, you only need to increase the total sample size.

Priors

Variable Definition
p.alpha this is actually log(alpha) in JAGS implementation ln.alpha ~ dnorm(p.alpha,tau_alpha) , but keeping spec file consistent with other implementations.
tau_alpha
p.beta this is actually the prior for capacity Smax, where S.max ~ dlnorm(p.beta, tau_beta)and beta <-1/S.max . The default is currently calculated as log( Max (Spn/sr.scale)) . A larger value for p.beta corresponds to a larger capacity, which translates into a smaller beta, which in turn means less of a density-dependent "penalty" on productivity in the Ricker form (see equations above)
tau_beta

| max.scalar | | |

We are currently exploring approaches for setting default priors based on the input values for a stock (e.g. observed range of spawners). For details, see discussion here.

TO INCORPORATE:

				mcmc.inits = "default",

				mcmc.priors = list(p.alpha = 0,tau_alpha = 0.0001, p.beta = 1 , tau_beta = 0.1,max.scalar = 3,
				                   shape.tau_R = 0.001,lambda_tau_R=0.01,shape.tauw = 0.01,lambda_tauw=0.001),

Clone this wiki locally