10-appendix-d.rmd

# Preparation of Data for Fitting Spawner-Recruit Models to Substocks of Kuskokwim River Chinook Salmon in Chapter \@ref(ch4) {#appendix-d}

## Overview of data needs

\noindent
All data for this analysis are available to the public, and came primarily from the Arctic-Yukon-Kuskokwim Database Management System (AYKDBMS)^[http://www.adfg.alaska.gov/CommFishR3/WebSite/AYKDBMSWebsite/Default.aspx] maintained by the Alaska Department of Fish and Game (ADF&G). Cases in which other data sources were necessary are highlighted in the description, _e_._g_., the telemetry data needed to perform the expansion of aerial survey counts described in Section \@ref(air-expansion) below.

\noindent
This analysis required three primary data sources:

(1)  Estimates of annual escapement to each of the substocks included.

(2)  Estimates of annual harvest. Linear regression models (Section \@ref(reg-methods)) required harvest apportioned to each substock, the state-space models (Section \@ref(ssm-model)) required only total aggregate harvest summed across all substocks included.

(3)  Estimates of annual age composition (_i_._e_., the fraction of the run each year made up of each age) for all substocks that have had it collected.

\noindent
Any of these data sources could have missing years. 

## Substock escapement {#air-expansion}

\noindent
Escapement count data for this analysis were informed predominately by the ADF&G Kuskokwim River salmon escapement monitoring program, the details of which have been most-recently documented in @head-smith-2018. The data set available spanned 20 different escapement monitoring projects (six weirs and 14 aerial surveys) and 42 calendar years from 1976 -- 2017, though monitoring projects differed in the date they were initialized (Figure \@ref(fig:obs-freq)). For substocks monitored _via_ weir, observed escapement for substock $j$ in year $t$ ($S_{obs,t,j}$) was taken to be the total estimated weir passage each year. Substocks monitored _via_ aerial survey needed special care, however. Surveys have been flown only once per year on a relatively small fraction of each tributary system (Figure \@ref(fig:ch4-map)), resulting in these data being indices of escapement rather than estimates of total escapement to the subdrainage. This analysis required estimates of total escapement to each substock however, because this would allow calculation of biological reference points that are expressed in terms of the scale of the population (_e_._g_., the spawner abundance that is expected to produce maximum recruitment;  $S_{\text{MAX},j}$), rather than as a rate (_i_._e_., $U_{\text{MSY},j}$).

The approach used to estimate total escapement from single-pass aerial surveys involved two main steps:

(1) Mapping the distribution of detected telemetry-tagged Chinook salmon against distribution of the aerial survey counts. This comparison allowed for a spatial expansion to estimate how many salmon would have been counted had the entire tributary been flown. Radio telemetry studies targeting Chinook salmon in the Kuskokwim River have been conducted intermittently since the early 2000s, and data were used from the years 2003 -- 2007 and 2015 -- 2016. These studies tagged fish migrating through the middle or lower river and (among other things) located them at the end of the season using aerial telemetry [@stuby-2007; @smith-liller-2017a; @smith-liller-2017b].  

(2) Obtaining and applying a temporal correction factor for the problem of counting a dynamic pool at one point in its trajectory. This correction factor was based on the relationship between paired weir and aerial counts on $n=3$ of the systems in the analysis.

### Spatial expansion {#spat-expansion}

\noindent
The core of the the spatial expansion estimator was the assumption:

\begin{equation}
  \frac{A_{f,t,i}}{T_{f,t,i}} = \frac{A_{u,t,i}}{T_{u,t,i}},
  (\#eq:air-expand1)
\end{equation}

\noindent
where the quantities $A$ and $T$ represent fish and tags, respectively, in flown ($A_f$ and $T_f$) and unflown ($A_u$ and $T_u$) reaches in year $t$ and for aerial survey monitoring project $i$. This assumption states that the ratio of actual spawners per one tagged spawner is the same between flown and unflown river sections at the time of the aerial index count and the aerial telemetry flights. Equation \@ref(eq:air-expand1) and can be rearranged as:

\begin{equation}
  A_{u,t,i} = A_{f,t,i} \frac{T_{u,t,j}}{T_{f,t,i}}.
  (\#eq:air-expand2)
\end{equation}

\noindent
If $T_{u,t,i}$ is further assumed to be a binomial random variable with time-constant success parameter $p_i$, then:

\begin{equation}
  T_{u,t,i} \sim \text{Binomial}(p_i,T_{u,t,i} + T_{f,t,i}).
  (\#eq:air-expand-binomial)
\end{equation}

\noindent
Here, $p_i$ represents the probability that a tagged fish in the spawning tributary monitored by project $i$ was outside of the survey flight reach at the time of the aerial telemetry flight. When \@ref(eq:air-expand-binomial) is  rearranged to put $p_i$ on the odds scale, then:

\begin{equation}
  \psi_i=\frac{p_i}{1-p_i}.
  (\#eq:air-expand-odds)
\end{equation}

\noindent
The odds value $\psi_i$ can be substituted for the division term in \@ref(eq:air-expand2) which gives:

\begin{equation}
  A_{u,t,i} = A_{f,t,i} \psi_i.
  (\#eq:air-expand3)
\end{equation}

\noindent
To obtain the total number of fish that would have been counted had the entire subdrainage been flown ($\hat{A}_{t,i}$), the components can be summed:

\begin{equation}
  \hat{A}_{t,i} = A_{f,t,i} + A_{u,t,i}.
  (\#eq:air-expand4)
\end{equation}

\noindent
Substitution of \@ref(eq:air-expand3) into \@ref(eq:air-expand4) and factoring gives the estimator:

\begin{equation}
  \hat{A}_{t,i}=A_{f,t,i}(1 + \psi_i).
  (\#eq:air-expand-final)
\end{equation}

\noindent
The spatial expansion model was integrated with the temporal expansion model described below into a single model fitted in the Bayesian framework fitted using JAGS [@plummer-2017]. This allowed for seamless propagation of uncertainty (in $\psi_i$) from the expansion above to the next step: a temporal expansion.

### Temporal Expansion {#temp-expansion}

\noindent
A temporal expansion model was necessary to convert from the one-pass index scale to the substock total annual escapement scale. The temporal expansion employed here operated by first regressing $n = 16$ observations of paired weir count ($W_i$) and spatially expanded aerial counts ($\hat{A}_{i}$; given by \@ref(eq:air-expand-final)) on the same tributary systems ($n = 3$) in the same years:

\begin{equation}
  \begin{split}
    W_i = \beta_0 + \beta_1 \hat{A}_i + \varepsilon_i, \\
    \varepsilon_i \stackrel{\text{iid}}{\sim} \text{N}(0, \sigma_W^2) \\
  \end{split}
(\#eq:temp-expand1)
\end{equation}

The estimated coefficients $\hat{\beta}_0$ and $\hat{\beta}_1$ (Table \@ref(tab:temp-expand-table)) were then applied to tributary systems with an aerial count but not a weir count:

\begin{equation}
  S_{obs,t,j}=\hat{\beta}_0 + \hat{\beta}_1 \hat{A}_{t,j}
(\#eq:temp-expand2)
\end{equation}

\noindent
The fitted relationship is shown in Figure \@ref(fig:obs-correct). For substocks that had both weirs and aerial surveys, the weir count was used as $S_{obs,t,j}$ as opposed to using the expansion in \@ref(eq:temp-expand2) and the coefficient of variation (CV) representing observation uncertainty for the state-space models was set at 5%, which assumed annual escapement counts made at weirs are made with little measurement error. For substocks monitored solely _via_ aerial survey, the posterior mean value of $S_{obs,t,j}$ was used as the escapement count that year, and the posterior CV was calculated for use as the observation uncertainty passed to the state-space models. In some cases, multiple aerial survey projects were considered to monitor a distinct portion of a larger substock (_e_._g_., three aerial survey projects count the three forks of the Aniak system; Figure \@ref(fig:ch4-map), Table \@ref(tab:spat-expand-table)). For this reason, expansions were conducted for each aerial survey project independently and then summed to obtain the estimates for that substock. 

## Aggregate harvest {#harv-expansion}

\noindent
Harvest estimates for the Kuskokwim River were available at the drainage-wide scale only, and were obtained each year by subtracting the drainage-wide estimates of total run and escapement [@liller-etal-2018]. Because the escapement data used here did not encompass all the substocks within the Kuskokwim River system, it was necessary to remove some portion of the total harvest that was produced by stocks not included in this analysis. First, the observed exploitation rate of the drainage-wide Kuskokwim River Chinook salmon stock ($U_{obs,t}$) was calculated by dividing the total harvest by the total run each year. Then, the assumption was made that monitored and unmonitored substocks have received the same exploitation rates, in which case total harvest accounted for in this analysis harvest could be obtained as:

\begin{equation}
  H_{obs,t} = \frac{S_{obs,t} U_{obs,t}}{1-U_{obs,t}},
  (\#eq:H-obs)
\end{equation}

\noindent
which can be derived from the definition of the exploitation rate $\left(U = \frac{H}{S+H}\right)$. This step was embedded within the same Bayesian model that encompassed the spatial and temporal aerial survey expansions such that uncertainty in these steps could be propagated through the entire analysis. The posterior mean value of $H_{obs,t}$ was used as the observed total harvest data, and the posterior CV was retained for use as the observation error attributed to this data source.

Note that $S_{obs,t}$ and $H_{obs,t}$ do not have $j$ subscripts denoting particular substocks: this indicates that they are aggregate quantities summed across all substock components. In cases where substock-specific harvest was required (_i_._e_., in reconstructing the substock-specific brood tables for fitting regression relationships; Appendix \@ref(lm-btable)), $H_{obs,t,j}$ was obtained using \@ref(eq:H-obs) by substituting $S_{obs,t,j}$ in for $S_{obs,t}$.

## Age composition {#age-comp}

\noindent
Age composition data were necessary to reconstruct brood tables for age-structured salmon populations (see Appendix \@ref(lm-btable)). Age data used in this analysis came from the ADF&G standardized age, sex, and length sampling program operated at the weir projects. All sampled fish that were not aged successfully were discarded as were samples corresponding to the rare ages of 3 and 8 such that only fish successfully aged as between 4 and 7 were included. It is possible that older or younger fish may have the systematic tendency to return early or late in the run, and this could introduce biases if age sampling was not conducted proportionally to fish passage throughout the season. To adjust for this possibility, a weighted-average scheme was applied to obtain the age composition estimates for each substock and year with data.  Daily age samples were stratified into two-week strata and strata-specific proportions-at-age were calculated. These strata-specific age compositions were then averaged across strata within a year and stock weighted by the number of Chinook salmon estimated to have passed the weir in each stratum. The total number of fish successfully aged for each year and substock was retained for data-weighting purposes for the state-space models, which (unlike the regression-based approaches) internally reconstructed the brood tables.

## Brood table reconstruction {#lm-btable}

\noindent
An important consideration in the use of the regression-based method (Section \@ref(reg-methods)) is in how the  $\text{RPS}_{y,j}$ data are obtained for salmon stocks that return at more than one age (like Chinook salmon), especially given the assumption of no observation error. Only the states $S_{y,j}$ are ever directly observed (and generally with some level of error); $R_{y,j}$ is observed (for Chinook salmon) over four calendar years as not all fish mature and make the spawning migration at the same age. Thus, in order to completely observe one $\text{RPS}_{y,j}$ outcome, escapement must be monitored in year $y$ and escapement, harvest, and age composition must be monitored in the subsequent years $y+4$, $y+5$, $y+6$, and $y+7$. It is evident that missing one year of sampling (which is common; Figure \@ref(fig:obs-freq)) can lead to issues with this approach. Only completely observed $\text{RPS}_{y,j}$ data were used for this analysis of reconstructing brood tables for regression analysis, with the exception of missing age composition data. For substocks with no age composition data (_i_._e_., those monitored _via_ aerial survey), the average age composition each year across substocks that have data was used to reconstruct $\text{RPS}_{y,j}$, but was provided only for years with escapement sampling for substock $j$. Only substocks with $\ge3$ completely observed pairs of $\text{RPS}_{y,j}$ and $S_{y,j}$ were included for model fitting, given each fitted line was dictated by two parameters.

## Results

### Escapement {#esc-data-results}

\noindent
The escapement project that received the largest spatial expansion factor ($1 + \hat{\psi}_j$) was the Holitna River, with an expansion factor of 4.78 (4.04 -- 5.73; 95% equal-tailed credible limits, Table \@ref(tab:spat-expand-table)). Given the small length of surveyed stream relative to the size of the Holitna subdrainage (Figure \@ref(fig:ch4-map); substock #7), this large estimate makes intuitive sense. The aerial survey project that required the smallest spatial expansion was the Salmon Fork of the Aniak River (1.04; 1.01 -- 1.14; Table \@ref(tab:spat-expand-table)); note that this project captures nearly all of this tributary of the Aniak system (western-most fork of drainage #4; Figure \@ref(fig:ch4-map)). The average spatial expansion factor across all aerial survey projects was 1.78 (1.67 -- 1.93). 

In terms of the temporal expansion, the estimate for the primary expansion coefficient ($\hat{\beta}_1$) was 2.3 (1.76 -- 2.85; Table \@ref(tab:temp-expand-table)). This estimate indicates that spatially corrected aerial survey counts needed to be scaled up by a factor of 2.3 in order to be consistent with the total annual escapement counts made at weirs in years and subdrainages that had paired counts. This estimate makes intuitive sense given ADF&G has structured the timing of aerial survey flights to coincide with the peak of the escapement timing arrival curve [@head-smith-2018], which should occur approximately halfway through the run, indicating the aerial counts would need to be doubled to account for the second half of the run. The freely estimated intercept of the temporal expansion ($\hat{\beta}_0$) had a value of 1.9 (-60.71 -- 62.4), indicating that if no fish are counted by an aerial survey, very little escapement actually occurred (Table \@ref(tab:temp-expand-table)). For the substocks monitored _via_ aerial survey that required these expansions, the average annual observation CV was estimated to be 18% (it was set at 5% for weir-monitored substock escapement).

Based on the scale of the summed escapement estimates from this analysis relative to the drainage-wide scale [which includes monitored and unmonitored substocks; @liller-etal-2018], approximately half (57%; 47% -- 66%) of the Chinook salmon aggregate population is accounted for by the $n_j$ = 13 substocks included here (Figure \@ref(fig:obs-fraction)).

### Harvest {#harv-data-results}

\noindent
The average annual exploitation rate used in this analysis was the same as for the aggregate Chinook salmon population: 0.42 with a minimum and maximum value of 0.13 and 0.62 in 2017 and 1988, respectively. The average annual harvest attributable to substocks included in this analysis was nearly 50,000, which is approximately 57% of the historical average harvest from all substocks in the Kuskokwim River. The average annual observation CV for total harvest was estimated to be 8%. 

### Age composition {#age-data-results}

\noindent
Age composition data were available for six substocks included in this analysis (_i_._e_., those monitored _via_ weir projects: George, Kogrukluk, Kwethluk, Takotna, Tatlawiksuk, and Tuluksak; Appendix \@ref(age-comp)). The number of years observed per substock depended on the starting year of the project: each year with an escapement count from a weir had associated age data but not all projects have been operating since 1976. The average annual multinomial sample size ($ESS_{t,j}$) varied among substocks (range 117 -- 553) and was generally related to the average escapement counted at each project. Across all substocks, the average annual age composition was 27%, 38%, 33%, and 2% for ages 4, 5, 6, and 7, respectively, though there was a high degree of inter-annual variability.

\clearpage

\singlespacing
```{r spat-expand-table}
ests = read.csv("tab/Ch4/spatial_expansion_estimates.csv")

ests = apply(ests, 2, add_break)

stks = c(
  "\\textbf{Kisaralik}", 
  paste(c("\\textbf{Salmon (Aniak)}", "\\textbf{Aniak}", "\\textbf{Kipchuk}"), footnote_marker_alphabet(1, "latex"), sep = ""),
  "\\textbf{Holokuk}", "\\textbf{Oskawalik}", "\\textbf{Holitna}", 
  paste(c("\\textbf{Cheeneetnuk}", "\\textbf{Gagaryah}"), footnote_marker_alphabet(2, "latex"), sep = ""),
  paste(c("\\textbf{Salmon (Pitka Fork)}", "\\textbf{Bear}", "\\textbf{Upper Pitka Fork}"), footnote_marker_alphabet(3, "latex"), sep = "")
)

ests = cbind(stks, ests)

colnames(ests) = c("\\textbf{Aerial Survey}", "$\\boldsymbol{\\hat{p}_i}$", "$\\boldsymbol{1 + \\hat{\\psi}_i}$")

kable(ests, "latex", booktabs = T, align = "lcc", 
      caption = "The estimated spatial expansion factors for the various aerial survey projects described in Appendix \\ref{spat-expansion}. $\\hat{p}_i$ represents the average fraction of telemetry tags that were detected outside of index flight reaches, which was used as the basis for determining the multiplier ($1 + \\hat{\\psi}_i$) needed to correct the aerial count for not flying the entire subdrainage. In cases where multiple projects were flown to count fish within one substock (\\textit{e}.\\textit{g}., the Aniak, see Figure \\ref{fig:ch4-map} substock \\#4), the expanded project counts were summed to obtain an estimate for the total substock, as indicated by the footnotes.",
      escape = F, linesep = "") %>%
  # column_spec(1, bold = T) %>%
  footnote(
    alphabet_title = "Substocks assessed with multiple aerial survey projects",
    alphabet = c("Aniak Substock", "Swift Substock", "Pitka Substock"),
    escape = F)
```

\clearpage

```{r temp-expand-table}

ests = read.csv("tab/Ch4/temporal_expansion_estimates.csv")

ests = apply(ests, 2, add_break)

param = c("$\\hat{\\beta_0}$", "$\\hat{\\beta_1}$", "$\\hat{\\sigma_W}$")

ests = cbind(param, ests)

colnames(ests) = c("\\textbf{Parameter}", "\\textbf{Estimate}")
rownames(ests) = NULL
kable(ests, "latex", booktabs = T, align = "lcc", 
      caption = "The estimated temporal expansion parameters for converting spatially expanded aerial counts to estimates of subdrainage-wide escapement abundance each year.",
      escape = F)
```

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch4/obs-freq.jpg}
  \caption{The frequency of escapement sampling for each substock monitored in the Kuskokwim River. Black points indicate years that were sampled for substocks monitored with a weir and grey points indicate years sampled for substocks monitored with aerial surveys. The vertical black line shows a break where > 50\% of the years were monitored for a stock.}
  \label{fig:obs-freq}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch4/obs-correct.png}
  \caption{The relationship between spatially expanded aerial survey estimates and weir counts during the same years and substocks as described by (\ref{eq:temp-expand1}). Notice the uncertainty expressed in the predictor variable; this was included in the analysis by incorporating both the spatial (Appendix \ref{spat-expansion}) and temporal (Appendix \ref{temp-expansion}) expansions in a single model fitted using Bayesian methods.}
  \label{fig:obs-correct}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch4/obs-fraction.png}
  \caption{Estimated Chinook salmon escapement for substocks within the Kuskokwim River drainage. ``Drainage-wide'' refers to the aggregate population estimates provided by a maximum likelihood run reconstruction model. ``This analysis'' refers to the estimated portion of the aggregate escapement included in this analysis (not all tributaries that produce Chinook salmon in the Kuskokwim River have been monitored; Figure \ref{fig:ch4-map}).}
  \label{fig:obs-fraction}
\end{figure}

\doublespacing