04-methods.Rmd

# METHODS  

We compiled available data on Fraser Pink Salmon spawner abundance and catch, then developed and fit a state-space spawner-recruitment model to these data to describe stock dynamics and population characteristics. 
We then derived estimates of biological reference points to assess stock status. 
Lastly, we developed a closed-loop simulation model conditioned on recent estimates of productivity to quantify future expected biological and fishery performance of the current HCR, an alternative HCR and a no fishing scenario for the stock. 
Each of these steps is described in detail below. 

## DATA SOURCES  

Spawner and catch data from 1959 to present were provided by the PSC [@pacificsalmoncommissionFraserRiverPanel2024]. 
Estimates of spawner abundance have been derived over this time period using a variety of approaches ranging from mark-recapture based methods in individual spawning tributaries or the Fraser mainstem to sonar based enumeration in the lower Fraser River in more recent years (Table \@ref(tab:tab-spawner-est-methods)). 
It should be noted that no calibration was done when switching between spawner abundance estimation methods, but methods prior to the extensive data review by @andrewReviewAssessmentAdult1987 have had corrections made, and more recent estimation methods benefited from lessons in that review and methodological advances. 
These approaches have varied in their precision but, based on conversations with area staff and other analysts familiar with the data, were assumed to account for the vast majority of the spawning population in any given year, and to not be systematically biased.  

Commercial catch has typically been estimated by multiplying the total coastal Pink catch (Canada and US) by the estimated contribution of Fraser stock to the coastal catch, where the contribution of the Fraser stock was estimated based on run-reconstructions (1959-85) and genetic stock identification methodologies (1987-present). 
Methods used to estimate catch vary by fishery type (i.e., commercial, recreational, First Nations) and country.
The Canadian commercial catch is estimated using a sales-slip program that began in 1951, while US catch is estimated through mandatory catch reporting to state fisheries departments in Washington and Oregon.
Data prior to `r min(rec_prop$year)` were not available because commercial landings did not partition catch into specific stocks, which is why our time series begins in `r min(rec_prop$year)`.
Recreational catch is estimated through agency creel surveys (e.g., effort surveys of catch per unit effort paired with flights or other counts) in both countries. 
First Nations Economic Opportunity catch is reported using similar methods of estimating commercial catch, while the First Nations Food, Social, and Ceremonial catch (FSC) is estimated using methods that differ by fishery and location. The total return (or recruitment) in any given year was assumed to be the sum of catch and spawner abundance.  

See @grantFraserRiverPink2014 for a detailed overview of the Fraser River Pink data landscape including methodologies used to collect spawner abundance, catch, and biological data. 

## SPAWNER-RECRUIT MODEL  

We modeled the spawner-recruitment data in a state-space framework, following the approach described in @fleischmanAgestructuredStatespaceStock2013. 
State-space models allow for the separation of observation (e.g., sampling) error and true underlying process variation and have become increasingly common in ecological modelling [@auger-metheGuideStateSpace2021]. 
State-space spawner-recruitment models tend to generate less biased estimates of leading parameters (e.g., intrinsic productivity and density dependence) than traditional regression based approaches that do not separate observation error and process variation and hence can be vulnerable to errors-in-variables and time series biases [@adkisonReviewSalmonSpawnerRecruitment2021; @statonEvaluationMethodsSpawner2020; @suPerformanceBayesianStatespace2012].  


### Process model  

The process model is intended to represent the true population dynamics (i.e., free of measurement error).
This component of our state-space spawner-recruitment model specifies productivity and density-dependence.
Recruitment abundances of adult Pinks ($R_y$) in odd year $y$ were treated as unobserved states and modeled as a function of spawner abundance in year ($S_{y-1}$) assuming a @rickerStockRecruitment1954 spawner-recruitment relationship with serially auto-correlated log-normal process variation:  

\begin{equation}
\ln(R_y) = \ln(S_{y-2}) + \ln(\alpha) - \beta S_{y-2} + v_y
\label{eq:AR1-ricker}
\end{equation}

where $\alpha$ is productivity (intrinsic rate of growth), $\beta$ is the magnitude of within brood year density-dependent effects, and $v_y$ reflects inter-annual variation in survival from egg to adulthood, which we term "recruitment anomalies".
This variation was assumed to follow a lag-1 autoregressive process (AR1) over time:

\begin{equation}
\begin{aligned}
v_y &= \phi v_{y-2} + \varepsilon_y \\
\varepsilon_y &\sim \mathcal{N}(0, \sigma_R)
\end{aligned}
\label{eq:AR1}
\end{equation}

where $\phi$ is the correlation coefficient and $\varepsilon_y$ reflects the portion of the recruitment anomaly $v_y$ that is temporally independent (i.e., white noise).
The first year of recruitment was not linked to observations of spawner abundance in the spawner-recruitment relationship (Eqn. \@ref(eq:AR1-ricker)) and were modeled as random draws from a log-normal distribution with mean $\ln(R_0)$ and standard deviation $\sigma_{R}^2$.
Rather than estimating $\ln(R_0)$ as a free parameter as in @fleischmanAgestructuredStatespaceStock2013, we choose to follow @statonEvaluationMethodsSpawner2020 and inform its value using the expected recruitment under equilibrium unfished conditions $\ln(\alpha)/\beta$.

Catch in a given odd year ($C_y$) was modeled as the product of total run size and the catch rate ($U_y$) experienced that year:

\begin{equation}
 C_y = R_y U_y
\label{eq:harvest}
\end{equation}

and spawner abundance ($S_y$) was modeled as the portion of $R_y$ remaining after catch $C_y$:

\begin{equation}
S_y = R_y (1 - U_y)
\label{eq:get-S}
\end{equation}


### Observation model  

We assumed that observation error in spawning abundance varied among assessment regimes, $r$ (Table \@ref(tab:tab-spawner-est-methods)):   

\begin{equation}
S_y = S_{obs_y} + \sigma^2_{r,y}
\label{eq:get-S}
\end{equation}

and then directly accounted for this by assuming observed spawner abundance was log-normally distributed with the coefficient of variation (CV) converted to log-normal variance following [@forbes_statistical_2011]:

\begin{equation}
\sigma^2_{r,y} = \ln\left(\mathrm{CV}_{r,y}^2 + 1\right)
\label{eq:get-sigma}
\end{equation}

We assumed that catch had a 5% CV and so catch observations were also assumed to be log-normally distributed with the CV converted to log-normal variance as per Eqn. \@ref(eq:get-sigma), then substituting catch, $C$, for spawners, $S$, into Eqn. \@ref(eq:get-S) and dropping the regime script, $r$.     


### Model fitting and diagnostics  

We fit the spawner-recruitment model in a Bayesian estimation framework with Stan [@carpenter_stan_2017; @standevelopmentteamRstanInterfaceStan2023], which implements the No-U-Turn Hamiltonian Markov chain Monte Carlo (MCMC) algorithm [@hoffman2014] for Bayesian statistical inference to generate a joint posterior probability distribution of all unknowns in the model. 
We sampled from 4 chains with 2,000 iterations each and discarded the first half as warm-up. 
We assessed chain convergence visually via trace plots and by ensuring that $\hat{R}$ (potential scale reduction factor; @vehtari2021rank) was less than 1.01 and that the effective sample size was greater than 200, or 10% or the iterations. 
Posterior predictive checks were used to make sure that the model returned data similar to the data used to fit parameters.  

Priors were generally uninformative or weakly informative and are summarized in Table \@ref(tab:tab-priors). 
The $\beta$ prior was moderately informative with a mean and variance of 75% of the maximum observed spawners, which prevents the model from exploring unrealistic parameter spaces of carrying capacity for Pacific Salmon (D. Greenberg, DFO, Nanaimo, British Columbia, pers. comm.).


## BIOLOGICAL REFERENCE POINTS  

We calculated biological reference points for each MCMC sample to propagate uncertainty. The spawning abundance expected to maximize sustainable yield over the long-term under equilibrium conditions, $S_{MSY}$ was derived as:  

\begin{equation}
S_{MSY} = 1 - W(e^{1-ln(\alpha)})/\beta
\label{eq:get-Smsy}
\end{equation}

where $W$ is the Lambert function [@scheuerellExplicitSolutionCalculating2016], and $\alpha$ and $\beta$ are intrinsic productivity and the magnitude of within stock density dependence, respectively. 
We chose to apply this exact solution for $S_{MSY}$ instead of the commonly applied @hilborn1985simplified approximation because the approximation only holds for $0 <ln(\alpha) \leq3$ such that infrequent, but large, posterior samples of $\alpha$ can result in biased estimates of the posterior distribution of $S_{MSY}$. 
We used 80% of $S_{MSY}$ as the USR following @holtEvaluationBenchmarksConservation2009 and @dfoSustainableFisheriesFramework2022.  

The catch rate expected to lead to maximum sustainable yield, $U_{MSY}$ was used as the RR and derived according to the solution proposed by @scheuerellExplicitSolutionCalculating2016 as:  

\begin{equation}
U_{MSY} = 1 - W(e^{1-ln(\alpha)})
\label{eq:get-Umsy}
\end{equation}

and $S_{gen}$, the spawner abundance expected to result in the stock rebuilding to $S_{MSY}$ in one generation in the absence of fishing [@holtEvaluationBenchmarksConservation2009], which we considered the LRP, was solved numerically according to: 

\begin{equation}
S_{MSY} = S_{gen}\alpha e^{-\beta S_{gen}}
\label{eq:get-Sgen}
\end{equation}

Equilibrium spawner abundance ($S_{eq}$), where recruitment exactly replaces spawners, was estimated as:  

\begin{equation}
S_{eq} = ln(\alpha)/\beta
\label{eq:get-Seq}
\end{equation}

## CLOSED-LOOP SIMULATION MODEL   

We developed a simple closed loop forward simulation, conditioned on our estimates of historical spawner abundance, and biological benchmarks illustrated in Figure  \@ref(fig:fig-schematic). 
We used this simulation to project the stock forward in time and evaluate the biological and fishery performance of the current, and an illustrative alternative HCR. 
Details on model components and calculation of performance are provided below.   


### Biological sub-model  

Because recruitment residuals tended to be negative in recent generations (Figure \@ref(fig:fig-rec-resid)), and reproductive potential has likely declined through time due to declining size (Figure \@ref(fig:fig-avg-mass)), we chose to refit a version of the model described in equation \@ref(eq:AR1) with time-varying intrinsic productivity that could then be used to condition the biological sub-model for the forward simulation. Specifically, we allowed the $\alpha$ parameter to evolve through time as a random walk, yielding annual estimates of productivity: 

\begin{equation}
\begin{aligned}
\alpha_y &= \alpha_{y-2} + \varepsilon_y \\
\varepsilon_y &\sim \mathcal{N}(0, \sigma_\alpha)
\end{aligned}
\label{eq:tv-alpha}
\end{equation}

and where recruitment anomalies were no longer modeled as being auto-correlated but all other parameters in equation \@ref(eq:AR1) otherwise remained the same.
We simulated future stock trajectories by starting with the most recent estimate (i.e., latent state) of spawners and the median estimate of productivity in the last 3 generations, then iterating the process model forward in time for five Pink Salmon generations (10 years) 
This was done 1000 times to ensure that uncertainty in the spawner-recruitment relationships was propagated by drawing for the joint posterior distributions of estimated parameters in each iteration of the simulation. 
 

### Fishery sub-model  

In each odd-year of the simulation, forecasted total returns of Pink Salmon were assumed to be estimated with error. 
This error was assumed to be lognormally distributed with a mean equal to the true run size and CV of `r round(for.error*100, 0)`% based on retrospective assessment of pre-season forecasts provided by the PSC for years 1987-2021. 
Forecasted returns were then used as an input into the HCR that specified target exploitation rate given the expected run-size. 
Outcome uncertainty (i..e., deviations from targeted catch) was then applied to calculate realized catch and spawning abundance 
We assumed this outcome uncertainty was normally distributed around the target catch with a CV of `r OU.CV*100`%.  

In addition to evaluating the current HCR we also considered an illustrative alternative and a no fishing scenario. 
The alternative HCR is fit to Fisheries and Oceans Canada's Precautionary Approach Framework [@dfoFisheryDecisionmakingFramework2009] (Figure \@ref(fig:fig-HCRs)). 
Under this alternative (“PA alternate”) HCR the lower operational control point (OCP) is set to our median estimate of $S_{gen}$ below which the target exploitation rate is zero, and an upper OCP is set to run-size associated with our median estimate of 80% $S_{MSY}$ where the maximum target exploitation rate is set to our median estimate of the RR or $U_{MSY}$. 
At run-sizes between the lower and upper OCPs the target exploitation rate was linearly interpolated. 
An unintended consequence of this alternative is that the associated target spawner abundance declines slightly as run-sizes increases toward the upper OCP which is undesirable and problematic from a practical management implementation perspective. 
To avoid this, slight variations on this type of HCR have been used for Fraser Sockeye [@pestalUpdatedMethodsAssessing2012].  


### Performance measures  

We quantified the expected performance of the HCRs against biological and fishery objectives and associated quantitative performance measures (Table \@ref(tab:tab-perf-metrics-descriptions)).
Biological objectives included minimizing the probability spawner abundances fall below the LRP ($S_{gen}$), and maximizing the probability the stock maintains spawner abundances above the USR (80% $S_{MSY}$) and is hence in a “healthy” or desirable state. 
The values used for these biological reference points were based on the joint posterior distributions of estimated parameters in each iteration of the simulation thereby ensuring that uncertainty in these reference points was explicitly accounted for in the performance measure calculations. 
The percentages reported are simply the percentage of simulation-years that fall above or below a reference point (e.g. if a single year of spawners within a simulation is above or below the reference point it is counted).  

Fishery objectives included maximizing average catch and inter-annual stability in catch, and maximizing the probability annual catch fall above a minimum catch index level which, for illustrative purposes, was chosen as the mean catch of the highest three catches since the year 2001 (i.e., an indicator of a "good" fishing year).  


### Robustness test  

We used a robustness test to evaluate the sensitivity of HCR performance to a potential situation where future intrinsic productivity of the Fraser Pink stock dramatically further declined due to, for example, large changes in marine or freshwater survival. To implement this we simply conditioned our biological sub-model using draws from the joint posterior distribution of Ricker parameters (i.e., $\alpha$, $\beta$, $\sigma$) associated with the lower tenth percentile of the median posterior distribution of the last 3 generations of the productivity parameter ($\alpha$), while draws from the starting state (i.e., spawners in 2023) and benchmarks were sampled from the full posterior distribution.