Conditioning model from CWT and escapement
2026-04-26
Source:vignettes/conditioning.Rmd
conditioning.RmdIntroduction
Typical run reconstruction uses a historical time series of spawner and recruit data to estimate the productivity of salmon populations. For many small populations of interest, escapement time series may be available but catch may not be identifiable to the local scale, for example, if the population is a conservation unit (CU) as part of a larger stock complex and catch composition is not identified to individual CUs.
Walters and Korman (2024) demonstrated an approach for reconstruction if Coded Wire Tag (CWT) data from an indicator hatchery are assumed to be representative to life cycle parameters of a natural population.
The model consists of two components. First, CWT data informs natural mortality, maturation, and exploitation rate in the marine environment. Second, these parameters are then applied to the population of interest and informs the size and productivity of the population that has a time series of total escapement. The number of hatchery releases in the system informs hatchery production, and the difference in total escapement and hatchery production informs natural production. Both steps are accomplished within a single model fit, which can account for uncertainty among the various data components, and posterior distributions of parameters are obtained by MCMC.
We utilize this approach as a conditioning model to inform stochastic parameters for projections in salmonMSE, although use of this model is not necessary to set up an operating model.
Model fitting
Model fitting is performed in RTMB with fit_CM().
The posterior can then be sampled with sample_CM()
through rstan.
A markdown report is available with report_CM()
A subset of posterior MCMC draws of parameters to reconstruct the
historical population can be imported into the operating model with
CM2SOM(). Further modification of the operating model and
additional settings can then be added to run the projection.
Variable definitions
| Name | Definition |
|---|---|
| \(C\) | Catch |
| \(\textrm{CWT}\) | Coded wire tag |
| \(E\) | Escapement |
| \(F\) | Instantaneous fishing mortality |
| \(\textrm{HO}\) | Hatchery origin |
| \(\textrm{HOS}\) | Hatchery origin spawner |
| \(M\) | Instantaneous natural mortality |
| \(N^\textrm{juv}\) | Juvenile abundance |
| \(\textrm{NO}\) | Natural origin |
| \(\textrm{NOS}\) | Natural origin spawner |
| \(\textrm{PT}\) | Preterminal fishery |
| \(R\) | Recruitment |
| \(\textrm{T}\) | Terminal fishery |
| \(m\) | Proportion mature |
| \(p^\textrm{female}\) | Proportion female |
| \(p^\textrm{spawn}\) | Proportion spawning |
| \(s^\textrm{enroute}\) | En-route survival (escapement to spawning grounds) |
| \(v\) | Fishery vulnerability |
Life cycle parameters from CWT
Life cycle parameters are informed by CWT. Mortality rates are parameterized in instantaneous units, which can be converted to survival terms.
Fishing mortality \(F\) is separated into a preterminal (PT) component which acts on juvenile fish and a terminal (T) component on mature fish. Separable effects are modeled where the fishing mortality is year-specific and modified by age class by a vulnerability term \(v\).
Survival of juvenile CWT to the next age class \(a\) at the beginning of year \(y\) is calculated after exploitation, maturation \(m\), and natural mortality \(M\):
\[ N^\textrm{juv,CWT}_{y+1,a+1} = N^\textrm{juv,CWT}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)(1 - m_{y,a})\exp(-M_{y,a}) \]
An additional mortality term, associated with release mortality, can be applied for survival in the first year of life shortly after release:
\[ N^\textrm{juv,CWT}_{y,a=1} = N^\textrm{rel,CWT}_y s^\textrm{rel} \]
where \(s^\textrm{rel}\) is the survival after release into the wild.
The CWT return \(R\) is the fraction of maturing juveniles after preterminal exploitation
\[ R^\textrm{CWT}_{y,a} = N^\textrm{juv,CWT}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)m_{y,a} \]
The escapement is the return that survive terminal exploitation
\[ E^\textrm{CWT}_{y,a} = R^\textrm{CWT}_{y,a}\exp(-v^\textrm{T}_aF^\textrm{T}_y) \]
The catch \(C\) is
\[\begin{align} C^\textrm{CWT,PT}_{y,a} &= N^\textrm{juv,CWT}_{y,a}(1 - \exp(-v^\textrm{PT}_aF^\textrm{PT}_y))\\ C^\textrm{CWT,T}_{y,a} &= R^\textrm{CWT}_{y,a}(1 - \exp(-v^\textrm{T}_aF^\textrm{T}_y)) \end{align}\]
Population reconstruction
The model does separate accounting of natural-origin (\(\textrm{NO}\)) and hatchery-origin (\(\textrm{HO}\)) fish in the population of interest.
Juvenile life stage
The age-1 HO fish is obtained from the number of releases and an assumption about survival after release: \(N^\textrm{juv,HO}_{y,a=1} = N^\textrm{rel}_y s^\textrm{rel}\).
The age-1 NO fish is equal to the smolt production, usually considered to be the juveniles migrating out of the freshwater environment: \(N^\textrm{juv,NO}_{y,a=1} = \textrm{Smolt}_y\)
In the first year in the marine life stage, both HO and NO fish experience the natural mortality rate associated with \(M\). Additional mortality can be applied to HO fish associated with release mortality (not experienced by NO fish). In effect, since the CWT data inform hatchery survival in the first year, the release mortality assumption effectively allows the analyst to attribute a smaller proportion of the mortality experienced by hatchery fish to natural origin fish.
The abundance of juvenile fish \(N\), uses the same preterminal fishing mortality, maturity, and natural mortality paramaters estimated from the CWT: \[\begin{align} N^\textrm{juv,NO}_{y+1,a+1} &= N^\textrm{juv,NO}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)(1 - m_{y,a})\exp(-M_{y,a})\\ N^\textrm{juv,HO}_{y+1,a+1} &= N^\textrm{juv,HO}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)(1 - m_{y,a})\exp(-M_{y,a}) \end{align}\]
Adult life stage
The return \(R\), and escapement \(E\) is calculated from the terminal fishing mortality and maturity parameters estimated from the CWT:
\[\begin{align} R^\textrm{NO}_{y,a} &= N^\textrm{juv,NO}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)m_{y,a}\\ R^\textrm{HO}_{y,a} &= N^\textrm{juv,HO}_{y,a}\exp(-v^\textrm{PT}_aF^\textrm{PT}_y)m_{y,a} \end{align}\]
\[\begin{align} E^\textrm{NO}_{y,a} &= R^\textrm{NO}_{y,a}\exp(-v^\textrm{T}_aF^\textrm{T}_y)\\ E^\textrm{HO}_{y,a} &= R^\textrm{HO}_{y,a}\exp(-v^\textrm{T}_aF^\textrm{T}_y) \end{align}\]
The number of spawners is calculated from the escapement, the en-route survival, and the proportion of escapement that spawn. This proportion can be calculated from the ratio of broodtake to observed escapement (the model is parameterized in this manner for numerical stability and it is assumed broodtake is non-selective with respect to origin):
\[\begin{align} \textrm{NOS}_{y,a} &= E^\textrm{NO}_{y,a} s^\textrm{enroute} p^\textrm{spawn}_y\\ \textrm{HOS}_{y,a} &= E^\textrm{HO}_{y,a} s^\textrm{enroute} p^\textrm{spawn}_y \end{align}\]
Egg and juvenile production
The egg production is calculated from the proportion females and fecundity \(f\) at age:
\[ \textrm{Egg}_y = p^\textrm{female}\sum_a f_a (\textrm{NOS}_{y,a} + \gamma \times \textrm{HOS}_{y,a}) \]
The smolt production is calculated from a Ricker stock-recruit function:
\[ \textrm{Smolt}_{y+1} = \alpha\times\textrm{Egg}_y \times \exp(-\beta \times \textrm{Egg}_y) \exp(-\delta_y) \]
where \(\delta_y\) is an annual deviation from the density-dependent function, expressed as an instantaneous mortality rate. Here, smolts can be considered to be outmigrating juveniles. Thus, the egg-smolt survival occurs over freshwater life stages.
Ricker stock-recruit function
Walters and Korman (2024) show that the realized egg-smolt mortality rate can be defined as
\[ M_y^\textrm{egg} = M^\textrm{egg}_\textrm{min} + \beta \times \textrm{Egg}_y + \delta_y \]
where \(M^\textrm{egg}_\textrm{min} = -\log(\phi^{-1}) - \log(\kappa)\) is the density-independent egg mortality rate when the egg production approaches zero, \(\beta \times \textrm{Egg}_y\) is the density-dependent mortality term, and \(\delta_y\) is a year-specific deviation.
The Ricker \(\alpha\) parameter is \(\alpha = \kappa/\phi\) in units of juvenile-per-egg, where \(\kappa\) is productivity (return/spawner) and \(\phi\) is the egg-per-juvenile at replacement (corresponding to the one-to-one return per spawner line), with \(\phi = \sum_a \ell_a m^\textrm{base}_a f_a p^\textrm{female}\).
\(\ell_a\) is the juvenile survival at age:
\[ \ell_a = \begin{cases} 1 &, a = 1\\ \ell_{a-1} \exp(-M^\textrm{base}_{a-1}) (1 - m^\textrm{base}_{a-1}) &, a = 2, \ldots, A \end{cases} \]
The \(\beta\) parameter is \(\beta = \log(\kappa)/(\phi \times J^\textrm{rep})\) is in units of reciprocal eggs, where \(J^\textrm{rep}\) is the juvenile smolt production at replacement. The equilibrium spawners at replacement is \(S^\textrm{rep}\) and \(J^\textrm{rep} = S^\textrm{rep}/\tau\), where \(\tau = \sum_a \ell_a m^\textrm{base}_a p^\textrm{female}\) is the equilibrium female spawners per juvenile.
Realized productivity
Note that the replacement values described above correspond to initial maturity and natural mortality values. The realized replacement, along with the realized productivity and spawners at replacement, varies if either maturity and natural mortality are varying in the historical time series because these parameters change the egg production over the lifetime of a juvenile.
To obtain the realized productivity in a particular year, we can calculate the egg-per-juvenile \(\phi_y = \sum_a \ell_{y,a} m_{y,a} f_a p^\textrm{female}\) and spawner-per-juvenile \(\tau_y = \sum_a \ell_{y,a} m_{y,a} p^\textrm{female}\) from the historical maturity and mortality, where
\[ \ell_{y,a} = \begin{cases} 1 &, a = 1\\ \ell_{y,a-1} \exp(-M_{y,a-1}) (1 - m_{y,a-1}) &, a = 2, \ldots, A \end{cases} \]
The realized productivity is \(\kappa_y = \alpha \phi_y\) with corresponding \(S^\textrm{rep}_y = \log(\kappa_y)\times\tau_y/(\beta\phi_y)\) and \(E^\textrm{rep}_y = S^\textrm{rep}_y \phi_y / \tau_y\) (\(E^\textrm{rep}_y\) is the egg production at replacement).
The realized productivity can be less than one (and \(S^\textrm{rep}_y < 0\)) in an individual year. If juvenile marine mortality is high or maturity is late, then there are too few returns to replace the population.
Initial abundance
As a forward-projecting model, an assumption about the equilibrium age distribution in the first year of the model is required. This age distribution is calculated from the initial fishing mortality rate \(F^\textrm{init}\) of preterminal and terminal fisheries.
The initial equilibrium juvenile survival at age is:
\[ \ell^\textrm{Finit}_a = \begin{cases} 1 &, a = 1\\ \ell_{a-1} \exp(-v^\textrm{PT}_{a-1}F^\textrm{init,PT})\exp(-M^\textrm{base}_{a-1}) (1 - m^\textrm{base}_{a-1}) &, a = 2, \ldots, A \end{cases} \]
The juvenile abundance is back-calculated from an initial equilibrium spawners \(S^\textrm{init}\) and pHOS assumptions:
\[\begin{align} N^\textrm{juv,NO}_{y=1,a} &= \dfrac{S^\textrm{init}}{\tau^\textrm{init}} \times (1 - p\textrm{HOS}^\textrm{init}) \times \ell^\textrm{Finit}_a\\ N^\textrm{juv,HO}_{y=1,a} &= \dfrac{S^\textrm{init}}{\tau^\textrm{init}} \times p\textrm{HOS}^\textrm{init} \times \ell^\textrm{Finit}_a \end{align}\]
where spawner-per-juvenile \(\tau^\textrm{init} = \sum_a \ell^\textrm{Finit}_a \exp(-v^\textrm{PT}_aF^\textrm{init,PT})\exp(-v^\textrm{T}_aF^\textrm{init,T}) m^\textrm{base}_a p^\textrm{female}\).
Parameter estimation and priors
Fishing mortality
Year-specific fishing mortality is parameterized as
\[\begin{align} F^\textrm{PT}_y &= \exp(a^\textrm{PT}) F^\textrm{trend,PT}_y \exp(\omega^\textrm{FPT}_y)\\ F^\textrm{T}_y &= \exp(a^\textrm{T}) F^\textrm{trend,T}_y \exp(\omega^\textrm{FT}_y)\\ \end{align}\]
where \(F^\textrm{trend,PT}\) is a time series of relative exploitation provided by the analyst. The model estimates a scaling coefficient \(a\) and annual deviations \(\omega\) to estimate fishing mortality.
The prior for the annual deviations is
\[\begin{align} \omega^\textrm{FPT}_y &\sim N(0, \sigma_\textrm{FPT}^2)\\ \omega^\textrm{FT}_y &\sim N(0, \sigma_\textrm{FT}^2) \end{align}\]
with hyperpriors for \(\sigma_\textrm{FPT}^2\) and \(\sigma_\textrm{FT}^2\):
\[\begin{align} \sigma_\textrm{FPT} &\sim \textrm{Gamma}(2, 5)\\ \sigma_\textrm{FT} &\sim \textrm{Gamma}(2, 5) \end{align}\]
Vulnerability
Vulnerability at age are independent terms estimated in logit space:
\[\begin{align} \textrm{logit}(v^\textrm{PT}_a) &\sim N(\textrm{logit}(\mu^{vPT}_a), 1.6^2)\\ \textrm{logit}(v^\textrm{T}_a) &\sim N(\textrm{logit}(\mu^{vT}_a), 1.6^2) \end{align}\]
where \(\mu\) is the prior mean specified by the analyst. Vulnerability is fixed at zero and one for the age 1 and maximum age (\(A\)), respectively (\(v_1 = 0\) and \(v_A = 1\)).
A useful prior may be with \(\mu = 0.5\). When transformed to normal space, the prior density is relatively uniform between 0-1 with low density at the bounds:
x <- seq(-5, 5, 0.1)
f_x <- dnorm(x, 0, 1.6)
y <- plogis(x)
g_y <- f_x /(y * (1 - y)) # Prior density with Jacobian transformation
par(mfcol = c(1, 2), mar = c(5, 4, 1, 1))
plot(y, g_y, type = 'l', xlab = expression(v[a]), ylab = "Prior density")
plot(y, log(g_y), type = 'l', xlab = expression(v[a]), ylab = "log prior density")
Maturity
Maturity is estimated in logit space as deviations from base parameters provided by the analyst. The prior density function is Gaussian with separate standard deviation \(\sigma^m_a\) by age:
\[ \textrm{logit}(m_{y,a}) \sim N(\textrm{logit}(m_a^\textrm{base}), [\sigma^m_a]^2) \]
with hyperprior \(\sigma^m_a \sim \textrm{Gamma}(2, 5)\) (mode of 0.2, mean of 0.4, and standard deviation of 0.28).
x <- seq(0, 1.5, 0.01)
f_x <- dgamma(x, 2, scale = 1/5)
#x[which.max(f_x)] # Prior mode at 0.2
par(mfcol = c(1, 2), mar = c(5, 4, 1, 1), oma = c(0, 0, 2, 0))
plot(x, f_x, type = "l", ylab = "Prior density")
plot(x, log(f_x), type = "l", ylab = "log prior density")
mtext("Gamma(2, 5) prior distribution", outer = TRUE, line = 1)
Natural mortality
Natural mortality is parameterized as
\[ M_{y,a} = \begin{cases} M^\textrm{base}_a + \sum_iX_{i,y}\beta_i+M^\textrm{add} + \varepsilon_y &, a = 1\\ M^\textrm{base}_a + \sum_jX_{j,y}\beta_j &, a = 2, \ldots, A-1 \end{cases} \]
From base values provided by the analyst, year-specific mortality rates can be estimated from linear combination of environmental covariates \(X\) and estimated coefficients \(\beta\). Separate covariates are used for age-1 and age-2+ fish.
For age-1, an additional scalar \(M^\textrm{add}\) and annual deviations \(\varepsilon_y\) can be estimated from the base parameter.
Gaussion priors are used for \(\varepsilon_y\):
\[ \varepsilon_y \sim N(0, [\sigma^M]^2) \]
with hyperprior \(\sigma^M \sim \textrm{Gamma}(2, 5)\).
Likelihoods
Total escapement
The likelihood of the total escapement for the population of interest uses a lognormal distribution:
\[ \log(E_y) \sim N\left(\log\left(\sum_a(\hat{E}^\textrm{NO}_{y,a} + \hat{E}^\textrm{HO}_{y,a})\right), [\sigma^E]^2\right) \]
pHOS
The model can also be fitted to observations of pHOS (proportion hatchery origin spawners by brood year \(t\)) for the population of interest. A logistic normal distribution is used:
\[ \log\left(\dfrac{p\textrm{HOS}_t}{1 - p\textrm{HOS}_t}\right) \sim N\left(\log\left(\dfrac{\widehat{p\textrm{HOS}}_t}{1 - \widehat{p\textrm{HOS}}_t}\right), [\sigma^{p\textrm{HOS}}]^2\right) \]
where the predicted value is
\[ p\textrm{HOS}_t = \dfrac{\sum_a\textrm{HOS}_{t+a-1}\times f_{t+a-1}}{\sum_a([\textrm{NOS}_{t+a-1} + \textrm{HOS}_{t+a-1})\times f_{t+a-1}]} \]
CWT catch and escapement at age
The likelihood of the CWT at age uses a Poisson distribution:
\[\begin{align} C^\textrm{CWT,PT}_{y,a} &\sim \textrm{Poisson}\left(\dfrac{1}{\lambda}\times\hat{C}^\textrm{CWT,PT}_{y,a}\right)\\ C^\textrm{CWT,T}_{y,a} &\sim \textrm{Poisson}\left(\dfrac{1}{\lambda}\times\hat{C}^\textrm{CWT,T}_{y,a}\right)\\ E^\textrm{CWT}_{y,a} &\sim \textrm{Poisson}\left(\dfrac{1}{\lambda}\times\hat{E}^\textrm{CWT}_{y,a}\right) \end{align}\]
where the \(\wedge\) symbol denotes an estimate of the corresponding observed quantity and \(\lambda\) is the expansion factor.
The expansion factor is the scalar to ensure CWT data are representative of total recoveries (typically greater or equal to 1). For example, \(\lambda = 10\) may be appropriate if 10 percent of the fishery catch is sampled for coded wire tags. It is assumed that the CWT data are unexpanded numbers.
Likelihood re-weighting
CWT likelihood components can be reweighted relative to the total escapement with alternative expansion factors (higher expansion factors have lower weight with increasing variance for the Poisson distribution).
The data needed to be adjusted simultaneously to preserve the magnitude of predicted CWT catch and escapement between alternative fits. For example, if the true sampling rate is 0.1, then an initial fit may use \(\lambda = 10\), but the CWT can be up-weighted assuming a higher sampling rate of 0.5 (\(\lambda = 2\)). The analyst should multiply the CWT catch and escapement by 5. # Release strategies
The equations above assume that CWT come from a single release strategy, for example, fed fry, smolts, or yearlings (see James et al. 2021 for a review).
The conditioning model can accommodate multiple release strategies. The number of CWT releases are specified by brood year and release strategy, and the model is fitted to CWT catch by year, age, and release strategy.
Fishing and natural mortality rates are identical by year and age but maturity differs by release strategy.
The maturity is estimated by year and age for a particular release strategy indexed by \(r = 1\):
\[ \textrm{logit}(m_{y,a,r}) \sim N(\textrm{logit}(m_a^\textrm{base}), [\sigma^m_a]^2) \]
For all other release strategies \(r = 2,\dots,n_r\), logit deviations by age (constant with time) are estimated as fixed effects:
\[ \textrm{logit}(m_{y,a,r}) = \textrm{logit}(m_{y,a,1}) + \delta_{a,r} \]
The user then specifies the release strategy from which the maturity estimates then apply towards the population of interest.