02-rt-fcst.Rmd

# Development and Evaluation of a Migration Timing Forecast Model for Kuskokwim River Chinook Salmon {#ch2}

\singlespacing
\noindent
\footnotesize
\textit{Note}: The majority of this chapter has been published in \cite{staton-etal-2017-rtf}. Changes included herein since publication include (\textit{a}) an updated analysis with two additional years (2017 and 2018) (\textit{b}) some embellishment on the sliding climate window approach, (\textit{c}) several new tables and figures, and (\textit{d}) revision of all figures to reflect additional data. All qualitative inferences made herein are identical to those published in \cite{staton-etal-2017-rtf}.

\vspace{0.25cm}

\noindent
Additionally, a more comprehensive analysis of the utility of the run timing forecast model for interpreting in-season run abundance index data is presented in @staton-catalano-2019.
\singlespacing

\normalsize
\doublespacing

## Abstract {-}

\noindent
Annual variation in adult salmon migration timing makes the interpretation of in-season assessment data difficult, contributing to in-season uncertainty about run size. I developed and evaluated a run timing forecast model for the Kuskokwim River Chinook salmon stock, located in western Alaska, intended to aid in reducing this source of uncertainty. An objective and adaptive approach (using model-averaging and a sliding window algorithm to select predictive time periods, both calibrated annually) was adopted to deal with multi-dimensional selection of four climatic variables and was based entirely on predictive performance. Forecast cross-validation was used to evaluate the performance of three forecasting approaches: the null (_i_._e_., intercept only) model, the single model with the lowest mean absolute error, and a model-averaged forecast across 16 nested linear models. As of 2018, the null model had the lowest mean absolute error (2.7 days), although the model-averaged forecast performed as well or better than the null model in the majority of retrospective years and currently has a mean absolute error of 3.2 days. The model-averaged forecast had a consistent mean absolute error regardless of the type of year (_i_._e_., average or extreme early/late) the forecast was made for, which was not true of the null model. Additionally, the availability of the run timing forecast was not found to increase overall accuracy but did reduce uncertainty of in-season run assessments in relation to the null model, particularly early in the season.

## Introduction

\noindent
In-season management strategies for Pacific salmon _Oncorhynchus_ spp. fisheries rely heavily on indices of in-river abundance (_e_._g_., test fisheries, sonar counts, _etc_.) to inform harvest control rules that attempt to attain the balance of meeting pre-determined escapement objectives while allowing adequate opportunity for harvest [@catalano-jones-2014]. However, because indices of abundance are confounded by the phenology (_i_._e_., timing) of the migration, their interpretation is very difficult in-season. For example, smaller-than-average index values early in the season could be due to either a small run with average timing or a late large run, when interpreted in the context of historical years [@adkison-cunningham-2015]. This ultimately leads to great uncertainty about how much of the incoming run has passed, which is a key piece of information that dictates fishery harvest opportunities. There exists no information in the current year’s abundance index to inform the manager if (for example) 25% or 75% of the run has passed on any given day. Yet, depending which is true, the optimal management decision could be vastly different. Thus, in-season assessment typically involves some characterization of the variation in historical run timing to formulate a range of possible run size scenarios that could be representative of the current year’s run size. However, given the amount of variation in historical run timing, these scenarios are rarely informative during the majority of the migration, when key harvest decisions are being made because the run scenarios may span all possible run sizes. As a result, the pre-season run size forecast may remain the most precise source of information for much of the season. If it were possible to predict the timing of the incoming run (_e_._g_., earlier- or later-than-average) with some level of confidence, it could prove valuable for in-season assessment and decision-making by reducing uncertainty in run size predictions.

While previous research has uncovered several key physiological mechanisms that are involved with natal homing [@hasler-scholz-1983] and return migrations of adult salmon to freshwater environments [@cooke-etal-2008; @cooperman-etal-2010; @hinch-etal-2012], the exact physiological and behavioral responses of adult salmon to relatively small-scale environmental gradients within estuaries, which are likely the ultimate determinants of freshwater entry timing, are still poorly understood. Despite this uncertainty, several hypotheses have been put forth that are broadly consistent with the observed timing patterns of several species across a large geographic area (_i_._e_., western and southwestern Alaska). Two primary influences have been suggested: genetic [@anderson-beer-2009; @omalley-etal-2010; @quinn-etal-2000] and environmental [@hodgson-etal-2006; @keefer-etal-2008] mechanisms. Substantial evidence exists to suggest that both genetic and environmental controls are involved in determining migration timing, however it is broadly thought that genetic variation influences substock variation (_i_._e_., different tributary spawning groups within the same major river basin) and environmental variation influences the timing of the aggregate (_i_._e_., basin-wide) run [@anderson-beer-2009; @keefer-etal-2008]. This is consistent with the notion that genetically distinct components of the aggregate run behave differently as a result of their life history strategies and/or the characteristics of their specific spawning grounds [_e_._g_., substocks that must travel farther in-river to reach spawning grounds enter freshwater earlier; @clark-etal-2015; substocks that spawn in tributaries influenced by warmer lakes enable later spawning; @burger-etal-1985] but that certain environmental conditions act on the aggregate run to either hasten or delay freshwater entry in any given year. It has also been suggested that run size may have an influence on migration timing, although empirical support for this claim seems to be lacking. If there were indeed relationships between run timing and run size, these need to be quantified as certain combinations are particularly troublesome for managers [_e_._g_., small/early runs and large/late runs appear the same early in-season; @adkison-cunningham-2015].

At the aggregate population scale, which is the focus of this Chapter, it has been observed that migrations occurring in the spring and summer generally occur earlier in years with warmer spring temperatures [@mundy-evenson-2011; @hodgson-etal-2006]. @mundy-evenson-2011 suggested that this pattern may be explained by the stability of the estuarine water column where adult salmon stage in preparation for riverine entry (or alternatively, marine exit). High estuarine water column stability was hypothesized to impede riverine entry through two mechanisms: (1) by presenting an osmotic barrier between freshwater riverine discharge and the saline ocean water which prevents osmotically incompetent individuals from crossing and, (2) by preventing freshwater-competent individuals from receiving olfactory cues essential to the homeward migration. Thus, @mundy-evenson-2011 hypothesized that years in which the estuarine water column is stable over a longer period of time would be associated with later migration timing. Although water column stability is a difficult variable to measure over large spatial scales, several variables that are known to influence it are available at large scales _via_ remote sensing (_e_._g_., satellite measurements). Such variables are sea ice cover which prevents wind-driven mixing, associated local temperature-related variables like land-based air temperature or sea surface temperature (SST), and broader scale indicators such as the Pacific Decadal Oscillation (PDO), an index of temperature anomalies in the northern Pacific Ocean. Observational studies across the North American range of Chinook salmon _O. tshawytscha_ have found environmental-run timing correlations that are consistent with this hypothesis [@hodgson-etal-2006; @keefer-etal-2008; @mundy-evenson-2011]. Even if the water column stability hypothesis is incorrect, observed patterns suggest that environmental variables may be useful in forecasting run timing with some level of accuracy and certainty; ultimately, understanding of the physical, physiological, or ecological mechanisms mediating run timing variability is of lesser importance than making good harvest and population management decisions.

Thus, several efforts have been made at exploiting these environmental-run timing relationships to develop run timing forecast models for Pacific salmon migrations. @mundy-evenson-2011 developed a model for Yukon River Chinook salmon that used air temperature, SST, and sea ice cover to predict the day at which the 15^th^ and 50^th^ percentiles of the run would pass a test fishery index location. Their model fitted the observed data well (errors were nearly always within seven days, usually within three days), although out-of-sample predictive ability was not presented. @keefer-etal-2008 developed a similar framework for Columbia River spring run Chinook salmon and found run timing relationships with river discharge, river temperature, and ocean condition indices (_e_._g_., PDO). Their best model explained 49% of the variation in median run timing with variation in the environmental variables. @anderson-beer-2009 continued this work on the Columbia River spring Chinook stock, but added genetic components to their analysis based on the arrival timing of precocious males. Their findings revealed that both environmental variables and substock composition affected overall run timing of the aggregate spring Chinook salmon run in the Columbia River. These advancements have shown that relationships between migration timing and environmental variables exist and may have utility for use in forecasting applications.

The Kuskokwim River, located in western Alaska, is the second largest river system in the state and supports culturally and economically important Chinook salmon fisheries. Chinook salmon return beginning in late May and continue through early August, with the median date of passage occurring between June 14 and July 2. Fisheries within the region harvest salmon in-river during freshwater migrations using primarily drift gillnet gear. The Kuskokwim River salmon fishery has a distinct cultural importance: nearly all inhabitants are native Alaskans belonging to the Yup’ik group and take salmon annually in the summer months for subsistence purposes [@linderman-bergstrom-2009]. While commercial salmon fisheries operate within the river, these fishers often also participate in subsistence take and revenues from the sale of commercially harvested salmon often contribute directly to participation in subsistence activities [@wolfe-spaeder-2009]. To ensure long-term sustainable harvest, the Chinook salmon fishery is managed with a drainage-wide escapement goal derived from an age-structured state-space spawner-recruit analysis [@hamazaki-etal-2012; @staton-etal-2017-intseq]. To meet these pre-determined escapement goals, in-season management strategies implement time, gear, and area closures based on limited and imprecise information regarding annual run size. The distant locations of the majority of escapement assessment projects makes direct measurement of escapement performance unavailable until late in the season. Thus, the primary sources of run size assessment information are (1) a pre-season run size forecast range (obtained as the previous year’s run size estimate $\pm \sim 20\%$) and (2) an in-river drift gillnet test fishery operated in Bethel, AK which has been implemented using consistent methods since 1984 [described in @bue-lipka-2016]. The interpretation of this test fishery index suffers from the same issue of being confounded by run timing described earlier, making management decisions difficult. Without precise in-season indicators of run size, managers must often choose to either trust a pre-season run size forecast for the majority of the season or somehow place weights on the various run timing hypotheses when interpreting in-season data. Both options could lead to the wrong interpretation of the actual run size, which could have serious consequences for the management of the fishery in a given year (_i_._e_., the unwarranted opening or closing the fishery resulting in severe under- or over-escapement). The necessity of more accurate and precise in-season perceptions of run size is particularly evident in years with anticipated low runs, such as in recent years (_i_._e_., since 2010 for Chinook salmon in the Kuskokwim and many other locations state-wide), as this may allow managers to more effectively guard against over-exploitation while still allowing for limited harvest opportunities to support the cultural and subsistence needs of the region. 

This chapter presents an analysis that develops and evaluates the performance of a run timing forecast model for Kuskokwim River Chinook salmon. The objectives were to:

(1)  quantify historical run timing, 
(2)  develop a run timing forecast model using environmental variables selected based on out-of-sample predictive performance,
(3)  assess the utility of the forecasting model for improving predictions of end-of-season test fishery indices of run size,
(4)  determine if there is a relationship between run size and run timing for the Kuskokwim River Chinook salmon stock.

## Methods

### Estimates of migration timing

\noindent
In this analysis, the forecasted quantity that characterized migration timing was the day at which 50% of the run passed an index location (hereafter, $D_{50}$). To inform this quantity for each year in the analysis, daily catch per effort (CPE) data were used from the Bethel Test Fishery (BTF) operated by the Alaska Department of Fish and Game (ADF&G), which spans the years 1984 -- 2018. The raw data were daily CPE beginning on June 1 and ending August 24 each year. The cumulative sum of these daily CPE values within a year follows a sigmoidal pattern reflecting the shape of the incoming salmon run which is characterized by relatively few early migrants, a peak where the majority of the fish are running, and relatively few late migrants. To estimate the median day of passage as a continuous variable, a logistic model was fitted to the cumulative proportion of daily CPE of the form:

\begin{equation}
  p_{d,t}=\frac{1}{1 + e^{-h_t (d - D_{50,t})}},
  (\#eq:logistic)
\end{equation}

\noindent
where $p_{d,t}$ is the predicted cumulative proportion on day-of-the-year (DOY) $d$ in calendar year $t$, $h_t$ is the parameter that controls the steepness of the curve (_i_._e_., duration of the run), and $D_{50,t}$ is the day at which 50% of the total annual CPE was caught in year $t$. Annual estimates of $D_{50,t}$ and $h_t$ were obtained by fitting $p_{d,t}$ to observed daily cumulative proportion by minimizing the sum of squared deviations from the model prediction. Statistical uncertainty from model fitting for these parameter estimates was not considered further in the analysis as the uncertainty was negligible (coefficients of variation for $D_{50}$ were less than 1%). Further, the use of the BTF daily CPE values to infer the location and shape of year-specific logistic timing curves made the assumption that these data provided an accurate representation of daily run strength within a year (_i_._e_., that the influence of weather conditions or harvest on sampling was negligible).

### Environmental variables

\noindent
Environmental variables assessed for forecasting performance were chosen based on three criteria: (1)  previously established association with salmon run timing, (2)  availability for the Kuskokwim River during the years for which BTF index observations exist (1984 – 2018), and (3)  availability for use in a pre-season forecast model (_i_._e_., available no later than June 10 in the year for which the forecasted value would be used). Based on these criteria, four environmental variables were chosen for analysis: SST, percent sea ice cover (SIC), PDO, and land-based air temperature taken in Bethel, AK.

#### PDO data

\noindent
Data collected for the PDO variable came from one of several indices produced by the National Oceanic and Atmospheric Administration (NOAA) [@mantua-etal-1997]^[PDO data: http://research.jisao.washington.edu/pdo/PDO.latest.txt]. The index is produced by taking the first principal component of monthly SST anomalies in the northern Pacific Ocean, after removing any global trends due to any systematic change over time [@mantua-etal-1997]. Thus, for each year of the data set, a single monthly value was available for PDO. Previous studies have found PDO values prior to the initiation of the run have predictive value for Chinook salmon populations [@beer-2007; @keefer-etal-2008].

#### Bethel air temperature data

\noindent
Air temperature data for Bethel, AK were accessed from the Alaska Climate Research Center^[Alaska air temperature data: http://akclimate.org/acis_data]. These data were available as daily means for each day of each year in the 1984 -- 2018 data set. @mundy-evenson-2011 found relationships with land-based air temperature and Chinook salmon run timing, which is why it was included in this analysis.

#### SST and SIC

\noindent
SST and SIC data were accessed from the NOAA Optimum Interpolation SST V2 High Resolution Dataset [@reynolds-etal-2007]^[Global gridded SST and SIC:  http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html]. These data were available as daily means for any 0.25° by 0.25° latitude by longitude grid cell on the globe. To limit the search, only grid cells within Kuskokwim Bay were selected for analysis (Figure \@ref(fig:ch2-map)) as that is the area where Chinook salmon bound for the Kuskokwim River likely aggregate prior to riverine entry. The area with grid cells ranged from 58.5° N to 60° N by 164.25° W to 162° W, which resulted in a total of 54 0.25° latitude by 0.25° longitude grid cells. For SST, four grid cells fell partially over land (resulting in 50 grid cells with daily data) and for SIC, five grid cells were partially over land (49 grid cells with daily data). “Empty” grid cells were excluded and the remaining grid cells were used for prediction. Previous analyses have used a simple average over a wide spatial area [_e_._g_., @mundy-evenson-2011] to obtain a single value for SST or SIC each year. However, this is somewhat arbitrary and does not account for the possibility of certain areas having stronger timing signals than others or that the areas with stronger signals may change over time. Thus, the gridded spatial structure of these variables was retained and the treatment of this structure in the forecast analysis is discussed below in Section \@ref(rtf-models).

### Forecast model {#reg-models}

\noindent
To produce a forecast of run timing, relationships between historical observed pairs of the environmental variables each year and $D_{50,t}$ must be quantified. The normal linear regression framework was used to obtain these historical relationships:

\begin{equation}
  \begin{split}
    D_{50,t} = \beta_0 + \beta_j x_{t,j} +, \cdots ,+ \beta_n x_{t,J} + \varepsilon_t, \\
    \varepsilon_t \stackrel{\text{iid}}{\sim} \text{N}(0, \sigma^2) \\
  \end{split}
(\#eq:lin-reg)
\end{equation}

\noindent
where $D_{50,t}$ is the observed run timing value in year $t$, $x_{t,j}$ is the observed value of covariate $j$ (of which there are $J$ included in the model), $\beta_0$ and $\beta_j$ are coefficients linking the observed values of $D_{50,t}$ with $x_{t,j}$, $\varepsilon_t$ are random residual effects that explain deviations of observed $D_{50,t}$ from the fitted value and have constant variance equal to $\sigma^2$. No time series patterns were visually present in the residuals and ARIMA models suggested no serial autocorrelation, which supports the use of the temporal independence assumption. 

There are many such regression models that could be used to produce a run timing forecast (_i_._e_., $\hat{D}_{50,t+1}$). This is because: (1) there are four variables (PDO, air temperature, SST, and SIC) that could be included, (2) each variable is temporally structured, _i_._e_., there are daily or monthly values for each variable, and (3) two variables (SST and SIC) are spatially structured, _i_._e_., there are different values for each day and year for different areas of Kuskokwim Bay (Figure \@ref(fig:ch2-map)). Item (1) deals with the specific values of $j$ and $J$ whereas items (2) and (3) deal with what values of $x_{t}$ for a given variable $j$ should take on. 

### Selection of predictive time periods {#clim-windows}

\noindent
To fit the regression model given by \@ref(eq:lin-reg), a single value for each $x_{t,j}$ was required. The covariate data were temporally structured, however, indicating that some selection of which time periods to use to populate $x_{t,j}$ was needed. Oftentimes the average over an arbitrary time period, such as daily values in the month of February, is used based on _a priori_ assumptions of the behavior of important factors [@vandepol-etal-2016]. While this approach is simple to implement and explain, it is possible that a better (_i_._e_., reliably more accurate predictor) time window exists but was not considered. Furthermore, the importance of various time windows may change over time and the arbitrary selection of a single window does not allow for such changes to be detected. To avoid these issues, a rigorous temporal selection process, known as the sliding climate window algorithm [SCWA; @vandepol-etal-2016], was implemented to determine the best predictive time period for each variable considered in the forecast model. To find the most reliable temporal window for prediction, the SCWA evaluates all possible windows (subject to certain restrictions) over which to average for use as the predictor variable in the forecast model. Although @vandepol-etal-2016 developed, discussed, applied, and simulation-tested the SCWA, no mathematical representation of the algorithm was presented. To provide more clarity on this approach to temporal variable selection, the following section provides the details of the SCWA.

#### The SCWA {#scwa}

\noindent
A "window" in this context is hereafter defined as a block of consecutive days in some portion of the year with starting day-of-the-year (DOY) denoted by $D_F$ and ending day equal to $D_L$. The daily values within each evaluated window (including $D_F$ and $D_L$) were averaged for the $x_{t,j}$ value to be used in a linear regression framework. As input constraints, the SCWA used in this analysis required: (1) the start DOY of the first window to be evaluated ($D_0$), (2) the end DOY of the last window to be evaluated ($D_n$), and (3) the minimum window size of a candidate window ($\Delta_{D,min}$). The algorithm started with the earliest and smallest possible time window: $D_F = D_0 = 1$ through $D_L = D_0 + \Delta_{D,min} - 1 = 5$. The performance of this window when used to obtain $x_{t,j}$ was evaluated for predictive performance (see Section \@ref(fcst-cv) below) and the result was stored for comparison to other candidate windows. For the next window, $D_F$ would remain at $D_0$, but $D_L$ would be incremented by 1 day ($\ell=1$). Thus, the endpoints of all candidate windows with $D_F = D_0$ can be generalized as:

\begin{equation}
  [D_0, D_0 + \Delta_{D,min} - 1 + \ell],
(\#eq:scwa-1)
\end{equation}

\noindent
for each $\ell = 0, 1, \cdots, n - 1$, where $n = D_L - D_F + 1$. For all windows, including those with $D_F = D_0$, this generalizes to:

\begin{equation}
  [D_0 + f, D_0 + f + \Delta_{D,min} - 1 + \ell],
(\#eq:scwa-2)
\end{equation}

\noindent
for each $f = 0, 1, \cdots, n - \Delta_{D,min}$ and $\ell = 0, 1, \cdots, n - \Delta_{D,min} - f$. Windows with $f > n - \Delta_{D,min}$ would contain fewer than $\Delta_{D,min}$ days and are thus prohibited according to the constraints. After evaluating all windows, the single window with the best predictive performance was used to obtain the forecast predictor variable for that data source (_i_._e_., PDO _versus_ air temperature).

As an example, consider the following inputs: $D_0 = 1$ (_i_._e_., January 1), $D_n = 31$ (_i_._e_., January 31), and $\Delta_{D,min}$ = 5. The SCWA would start with the window encompassing January 1 -- January 5, then do January 1 -- 6, January 1 -- 7, _etc_., January 1 -- 31. Next, it would exclude January 1 from consideration and evaluate all windows starting with January 2, starting with January 2 -- January 6. When it completes the one window starting with January 27, it must stop because windows starting later than January 27 would result in windows shorter than five days.

The values of $D_0$ and $D_L$ for the four covariates are shown in Table \@ref(tab:scwa-dates-table). The setting for $\Delta_{D,min}$  for air temperature, SST, and SIC was five days. Note that because PDO was available in monthly values only, each month was treated as "day" in the algorithm described above and $\Delta_{D,min}$ was set to one. 

#### Forecast cross-validation {#fcst-cv}

\noindent
A metric was needed to objectively measure the performance of the many evaluated windows. A time series forecast cross-validation procedure was used, which is an out-of-sample technique for data that are collected through time [@arlot-celisse-2010]. The procedure operated by producing a forecasted value of $D_{50}$ for year $t+1$ trained based on all data $x_{t,j}$ available from years $1, \cdots, t$. It then continued for all $t = m, \cdots, n-1$, where $m$ is the minimum number of years necessary to fit the model (set at $m = 10$ in all cases) and $n$ is the number of years of available data. Then, absolute forecast error was calculated based on all forecasted years as $|D_{50,t+1} - \hat{D}_{50,t+1}|$, and yearly forecast errors were averaged to obtain mean absolute error ($\overline{\text{AE}}$) which was used as the metric of model performance in window selection. The window with the lowest $\overline{\text{AE}}$ was selected as the optimal window to average over for prediction. The forecasting cross-validation procedure was used as opposed to other out-of-sample validation procedures, such as $k$-fold or leave-one-out methods, because the data were collected through time and the forecast model would never need to predict (for example) year 2010 from years 1984 -- 2009 and 2011 -- 2018, but rather it would always need to predict year $t+1$ from all previously collected data.

When forecasting $D_{50,t+1}$ from training data from $1, \cdots,t$, a single optimal climate window was selected for each variable and that window was used to estimate coefficients based on training data and obtain the environmental variable value for prediction in year $t+1$ to forecast $D_{50,t+1}$. When a new year of data was added to the training data (such as in the retrospective forecast analysis; Section \@ref(retro)), the optimal window for each variable was re-assessed using the algorithm again. For PDO and Bethel air temperature (which had no spatial structure), the SCWA was used to select the range of monthly (for PDO) or daily (for Bethel air temperature) values to include in the predictive climate window for each year in the analysis. For SST and SIC, the SCWA was used on each grid cell separately. The result was 50 unique grid cell-specific windows for SST and 49 windows for SIC for each year of the analysis. The treatment of this spatial structure in the forecast analysis is discussed below in Section \@ref(rtf-models).

### Evaluated forecast models {#rtf-models}

\noindent
Linear regression \@ref(eq:lin-reg) was used to assess the forecast performance of each of the variables described above, both in isolation of and in combination with other variables. All possible subsets were evaluated (excluding interactive effects) for predictive ability through time, resulting in a total of 16 models ranging from the null ($\beta_0$ only; $J = 0$) model to the full model (all $J=4$ variables included as additive predictors). 

For the spatially explicit variables (_i_._e_., SST and SIC), a more complex treatment was required to prevent all grid cell values from being used as predictors in a single model. To handle the spatial structure, grid cell-specific regression models were fitted, then model-averaging [@burnham-anderson-2002] based on Akaike's Information Criterion [AIC, @akaike-1974] was used to obtain a single forecast $D_{50}$ for each year. Under this approach, each grid cell $g$ received an AIC score:

\begin{equation}
  \text{AIC}_{\text{c},g}=n \log{\left(\hat{\sigma}_g^2\right) + 2K + \frac{2K(K+1)}{n-K-1}},
  (\#eq:aicc)
\end{equation}

\noindent
where $n$ is the number of data points used in each model, $\hat{\sigma}_g$  is the estimate of the residual standard deviation under grid $g$, and $K$ is the number of model parameters. The corrected version of AIC (AIC~c~) is recommended in cases where the ratio of $n$ to $K$ is small [@burnham-anderson-2002]. Then, each grid cell received a $\Delta\text{AIC}_\text{c}$ score, representing its relative performance in comparison to the best grid cell:

\begin{equation}
  \Delta_g=\text{AIC}_{\text{c},g}-\text{AIC}_{\text{c},min},
  (\#eq:delta-aicc)
\end{equation}

\noindent
where $\text{AIC}_{\text{c},min}$ is the minimum AIC~c~ across all grid cells. Model (grid cell) weights were then calculated as:

\begin{equation}
  w_g=\frac{e^{-0.5\Delta_g}}{\sum_j^G e^{-0.5\Delta_j}},
(\#eq:aicc-weights)
\end{equation}

\noindent
where $G$ is the number of grid cells. Grid cell-averaged predictions were then obtained as:

\begin{equation}
  \hat{y}_{t+1}=\sum_g^G w_g \hat{y}_{g,t+1},
(\#eq:grid-avg-fcst)
\end{equation}

\noindent
where $\hat{y}_{g,t+1}$ is the forecasted value of $D_{50}$ for grid cell $g$.

### Forecast uncertainty

\noindent
In addition to forecast accuracy, forecast precision is also of great importance as it is an index of the confidence with which the forecast should be interpreted. For models that did not require AIC~c~ model-averaging across grid cells, the following estimator was used to produce a forecast standard error (SE):

\begin{equation}
  \text{SE}=\hat{\sigma} \sqrt{1 + \frac{1}{n} + \frac{(x-\bar{x})^2}{\sum_i^n(x_i-\bar{x})^2}},
(\#eq:se)
\end{equation}

\noindent
where $n$ is the number of years the model was fitted to, $x$ is the value of the predictor variable used for forecasting, and $\bar{x}$ is the mean of all predictor values excluding the new value used for forecasting. For models that used AIC~c~ model-averaging (_i_._e_., those including SST and SIC), the following estimator was used to produce prediction SE:

\begin{equation}
  \text{SE}=\sum_g^G w_g \sqrt{\text{SE}_g^2+(\hat{y}_{g,t+1}-\hat{y}_{t+1})^2},
(\#eq:mod-avg-se)
\end{equation}

\noindent
where $\text{SE}_g$ is the prediction SE from grid cell $g$ calculated using \@ref(eq:se). This estimator of unconditional sampling standard error accounts for uncertainty within each model and the uncertainty due to model selection [@burnham-anderson-2002]. Prediction intervals were calculated using the point estimate of prediction, the prediction SE, and appropriate quantiles from the corresponding $t$-distribution. 

### Forecast model selection {#model-selection}

\noindent
Given 16 forecast models, it is impossible to know which may produce the best forecast for the current year. Thus, three methods to obtain a forecast for $D_{50}$ were evaluated: (1) the null (_i_._e_., intercept only) model, (2) the single model with the lowest forecast cross-validation score as of the last year, and (3) model-averaging across the ensemble of 16 forecast models based on AIC~c~ scores. According to @burnham-anderson-2002, model-averaging should perform better than a single "best model" at prediction when there is a high degree of uncertainty about which model is best. This procedure was performed using \@ref(eq:aicc) -- \@ref(eq:mod-avg-se), by substituting the prediction, prediction SE, and $K$ for forecast model $i$, in place of grid $g$. Prediction intervals based on model-averaged predictions and prediction SE present somewhat of a problem when the different models contributing to the average contain differing degrees of freedom as it is unclear how many standard errors the prediction limits should lie from the mean prediction. Thus, the estimator suggested by @burnham-anderson-2002 of the “adjusted SE” (ASE) was used:

\begin{equation}
  ASE=\sum_i^{16} w_i \sqrt{\left(\frac{t_{df_i,1-\alpha/2}}{z_{1-\alpha/2}}\right)^2 \text{SE}_i^2+(\hat{y}_{i,t+1}-\hat{y}_{t+1})^2},
(\#eq:ase)
\end{equation}

\noindent
where $t_{df,i,1-\alpha/2}$ is the $1-\alpha/2$ quantile of the $t$-distribution with degrees of freedom equal to that of model $i$ and $z_{1-\alpha/2}$ is the corresponding quantile of the $z$-distribution (_i_._e_., standard normal). The confidence level $\alpha = 0.05$ was used in all cases.

### Retrospective forecast analysis {#retro}

\noindent
The analysis was conducted in a retrospective forecast framework starting in 1994, which allowed the models to be trained to 10 years of observed data before being used for forecasting. All data after 1994 were ignored, optimal windows were selected for each of the four variables (and all grids for SST and SIC), all 16 models were fitted, a $D_{50}$ forecast was made for 1995 using the three approaches described in Section \@ref(model-selection), and each was evaluated for predictive accuracy. This process was repeated annually until the present (_i_._e_., out-of-sample predictions made for 1995 -- 2018), which allowed for the calculation of $\overline{\text{AE}}$ through time as if the forecast model would have been available beginning in spring 1995. In addition to $\overline{\text{AE}}$, median absolute error ($\widetilde{\text{AE}}$) was calculated to validate prediction accuracy of estimates by ignoring the effect of outlying poor predictions.

### Value of forecast to run size assessments {#eos-preds}

\noindent
It is important to remember that the purpose of producing a run timing forecast is to aid in the interpretation of in-season indices of run size such as test fisheries. To evaluate the utility of having access to the run timing forecast model, the accuracy and precision of an imperfect abundance index for the Kuskokwim River were compared when informed using $D_{50}$ forecasts from the model-averaged and the null forecast models. The abundance index is denoted by $\text{EOS}_t$, and is the end-of-season cumulative CPE observed in the BTF in year $t$. Under the assumption of constant catchability, $\text{EOS}_t$ should be proportional to total abundance, with deviations introduced by sampling noise. In-season predictions of $\text{EOS}_t$ were made for each year $t$, model $i$, and day $d$ in the season with: 

\begin{equation}
  \widehat{\text{EOS}}_{d,t,i}=\frac{\text{CCPE}_{d,t}}{\hat{p}_{d,t,i}},
(\#eq:eos)
\end{equation}

\noindent
where $\text{CCPE}_{d,t}$ is the cumulative CPE caught at the BTF through day $d$ in forecasting year $t$ $\left(\text{CCPE}_d = \sum_{j=1}^{d} \text{CPE}_j \right)$, $\hat{p}_{d,t,i}$ is the predicted cumulative proportion of the run that had passed the BTF location on day $d$ in year $t$ from model $i$ (_i_._e_., model-averaged _versus_ null forecast model) obtained by inserting the forecasted value of $D_{50}$ into the logistic function \@ref(eq:logistic). Uncertainty in the run timing forecast model was propagated to $\text{EOS}_{d,t,i}$ using a parametric Monte Carlo procedure [@bolker-2008]. Independent normal random samples for $D_{50,t}$ and $h_t$ were sampled to obtain $p_{d,t,i}$ for use in \@ref(eq:eos). The distribution for generating samples of $D_{50}$ differed based on the mean and SE of the forecasted value of $D_{50}$ according to either the null or the model-averaged model. Estimates of $D_{50,t}$ and $h_t$ seem to be independent for the Kuskokwim River Chinook stock, so the method of drawing samples of $h_t$ was identical for the comparison and involved sampling from a normal distribution with mean and standard deviation equal to the estimated quantities from all years before year $t$ to ensure the consistency of out-of-sample predictions. A sufficiently large number of Monte Carlo samples were drawn (10,000) for each evaluated day and year. Prediction uncertainty was quantified using the coefficient of variation (calculated as median/sample standard deviation across all Monte Carlo samples for year $t$ and day $d$).  Accuracy was assessed using mean absolute percent error (MAPE) and mean percent error (MPE), where the point estimate used was the median of all Monte Carlo samples of $\text{EOS}_{d,t,i}$. Prediction performance measures were compared between the null and the model-averaged forecast model on June 15, June 30, July 15, and July 30 each year a forecast was available (1995 -- 2018), which represent the historical average 21%, 79%, 95%, and 99% points of the Chinook salmon run in the Kuskokwim River. Using the null model to obtain $\hat{p}_{d,t}$ is one of several methods that could be done to produce the predictions $\text{EOS}_{d,t}$ in the absence of an environmental forecast variable model for $D_{50}$. This method was used (as opposed to other methods, like simply taking the average $p_d$ across all years to populate \@ref(eq:eos) so the comparison involved the same assumption of symmetry in the sigmoidal timing curve, which has the ability to affect the accuracy of $\text{EOS}_{d,t,i}$ predictions.

### Investigation of a run timing _versus_ run size relationship

\noindent
To test the hypothesis that run timing is related to run size (_e_._g_., small runs are typically early, or _vice versa_), two models were investigated for their predictive performance using the forecast cross-validation criteria: the null model and a model that included run size as a predictive covariate in place of the environmental variables. Run size was obtained from a maximum likelihood run reconstruction model that compiles all assessment information (_i_._e_., 20 escapement count indices, harvest estimates, drainage-wide mark-recapture estimates, _etc_.) to estimate the run size that makes the collected data most likely to have been observed [@bue-etal-2012; @liller-etal-2018]. The forecast absolute errors in each year were then compared using a two-tailed paired $t$-test with $\alpha = 0.05$.

## Results

### Estimates of run timing

\noindent
There was a considerable amount of inter-annual variability in $D_{50}$, with a range of 17 days and sample standard deviation of 3.62 days over the 35 years with run timing data from the BTF (Figure \@ref(fig:p-ccpue); Table \@ref(tab:rt-ests-table)). Based on $D_{50}$ alone, the earliest run on record was in 1996 when $D_{50}$ occurred on DOY 166 (June 14) and the latest run was in 1985 when $D_{50}$ was attained on DOY 183 (July 2). The average $D_{50}$ was 173.75 across all years available (which rounded down is June 22 in a normal year and June 21 in a leap year).

The logistic curve fitted the daily cumulative CPE proportions well in all years of the BTF data set (Table \@ref(tab:rt-ests-table)), as indicated by an average residual standard error estimate of 0.022, with a maximum estimate of 0.038 in 1992. The majority (95%) of all residuals from all years fell between -0.056 and 0.044. Parameter estimates were quite precise, with $D_{50}$ having a smaller average coefficient of variation than $h$, (0.07% and 2.07%, respectively). Given this negligible degree of parameter uncertainty, it was ignored throughout the rest of the analysis.

### Variable-specific relationships

\noindent
Looking at each of the environmental variables in isolation of all others, it is clear that there is a distinct relationship between temperature-related environmental variables and Kuskokwim River Chinook salmon migration timing (Figure \@ref(fig:relationships)). For illustration purposes, the figures for the two gridded variables (SST and SIC) were produced by taking an average across all grid cells weighted by the AIC~c~ weight for each grid cell. Air temperature, PDO, and SST all had negative relationships with $D_{50}$, whereas SIC had a positive relationship (Table \@ref(tab:coefs-table)). All single variable relationships were significant at the 0.01 level, but the estimated residual standard error was approximately 3 days for each model (Table \@ref(tab:coefs-table)). $R^2$ values were generally low (range: 0.19 for air temperature and 0.33 for SIC; Table \@ref(tab:coefs-table)).

### Selected climate windows

\noindent
It was difficult to generalize on the climate windows selected for each variable based on forecast cross-validation performance, because the selected windows changed with each new year of data and SST and SIC had windows for each grid cell; however, some noteworthy patterns arose. First, the best window for PDO was consistently the value for the month of May for each year the forecasts were produced (not shown). Second, selected windows for air temperature fluctuated from year to year to some extent: all short and mid-May windows were selected, then shifted to long windows spanning February to late-May beginning in the early 2000s (Figure \@ref(fig:window-changes)_a_). Third, selected windows through time were substantially more variable for most grid cells for SST and SIC than air temperature, although many grid cells remained relatively constant or became more "focused" as more years of data were added (Figures \@ref(fig:window-changes)_b_,_c_). In general, chosen windows for SST began in early- to mid-May and ended in late-May (Figure \@ref(fig:window-changes)_b_) whereas windows starting in early-April and ending in mid- to late-April were predominately chosen for SIC (Figure \@ref(fig:window-changes)_c_). The selected climate windows in southern-most grid cells appeared more stable for SST (Figure \@ref(fig:window-changes) panels \textit{b}3 and \textit{b}4), whereas climate windows in northern grid cells appeared more stable for SIC (Figure \@ref(fig:window-changes) panels \textit{c}1 and \textit{c}2; stable in the sense that the optimal windows changed less as new years were added to the training data). Finally, the best predictive windows tended to become systematically earlier in the season as the retrospective analysis progressed.

### Forecast performance

\noindent
Of the three investigated forecast methods (null model, model with lowest forecast cross-validation error up to the forecasting year, and AIC~c~ model-averaging), the null model had the lowest $\overline{\text{AE}}$ from 1995 to 2018 (2.7 days; Figure \@ref(fig:forecasts)). AIC~c~ model-averaging performed the same as using the single model with the lowest cross-validation score (both had $\overline{\text{AE}}$ = 3.2 days; Figure \@ref(fig:forecasts)). However, these patterns were not consistent across the entire time series. For the period of 1996 to 2008, the model-averaged forecast had a lower $\overline{\text{AE}}$ than the null model, and for the period of 2009 to 2015, the model-averaged forecast had approximately the same or lower $\overline{\text{AE}}$ scores (Figure \@ref(fig:ae-changes)). It was due in a large part to an extreme error in 2016 (Figure \@ref(fig:forecasts)) that the model-averaged forecast had a higher $\overline{\text{AE}}$ than the null model. A similar case occurred in 2015, but not 2014 when the run was truly very early (Figure \@ref(fig:forecasts)). Expressing prediction error in terms of median absolute error ($\widetilde{\text{AE}}$) resulted in lower average errors (null = 2.1, single best = 3.1, and model-averaged = 2.4), indicating that extreme prediction errors (_i_._e_., outliers) influenced the value of $\overline{\text{AE}}$ for the model-averaged forecast and the null model the most -- $\overline{\text{AE}}$ and $\widetilde{\text{AE}}$ were nearly identical for the approach that used the single model with the lowest CV score at the time (Figure \@ref(fig:forecasts)). Additionally, by comparing the width of the prediction intervals in Figure \@ref(fig:forecasts) across forecasting approaches, it was clear that model-averaging substantially reduced prediction uncertainty (SE) in relation to the null and single best model approaches.

To compare performance in average _versus_ extreme years among forecasting approaches, $\overline{\text{AE}}$ was further calculated in a more specific way: based on how similar or dissimilar the included years were to the mean observed run timing across all years. As would be expected, the null model performed well when years with $D_{50}$ within only $\pm 1$ days of the average were in included in the calculation of $\overline{\text{AE}}$ (Figure \@ref(fig:mae-subsets)_a_), but its accuracy became increasingly worse as years with more extreme realized $D_{50}$ values were included in the calculation (increasing $x$-axis values in Figure \@ref(fig:mae-subsets)_a_). The two environmental variable forecast approaches (model-averaging or the single "best" model in each year) performed nearly the same across this continuum and neither $\overline{\text{AE}}$ score was sensitive to the overall similarity or dissimilarity the included years had with average run timing (Figure \@ref(fig:mae-subsets)_a_). The lower panel shows the relative frequency with which these various scenarios occurred, indicating how much information each scenario contributed to the overall $\overline{\text{AE}}$. On the other hand, the null model only performed as well as the model-averaged forecast and single best model when years with $D_{50} \pm 0.5$ days outside of the mean were considered (Figure \@ref(fig:mae-subsets)_b_). As only more extreme years were considered in $\overline{\text{AE}}$ (increasing $x$-axis values in Figure \@ref(fig:mae-subsets)_b_), the null model rapidly performed worse and the model-averaged forecast remained relatively insensitive to the degree of extremity of the $D_{50}$ value it was forecasting relative to the mean (Figure \@ref(fig:mae-subsets)_b_). 

### Value to in-season run size assessments

\noindent
When the model-averaged and null model forecasts for $D_{50,t}$ were used to retrospectively assess the potential for the run timing forecast model to improve in-season run assessment based on daily cumulative BTF CPE, it was evident that the range of possible $\widehat{\text{EOS}}_{d,t}$ was substantially smaller when the model-averaged forecast was used as opposed to the null forecast. This is evident by the average daily coefficients of variation in the predictions made for $\widehat{\text{EOS}}_{d,t}$ on the two evaluated dates in June: 71% _versus_ 133% on June 15 and 14% _versus_ 22% on June 30 (Table \@ref(tab:eos-table)). The reduction in uncertainty of $\widehat{\text{EOS}}_{d,t}$ predictions in the first evaluated day is of importance as it is the time where many harvest management decisions are being made. Conversely, negligible improvements in the accuracy of $\widehat{\text{EOS}}_{d,t}$ predictions were found for the model-averaged forecast model in comparison to the null model. On June 15, MAPE when using $p_{d,t}$ informed by the model-averaged $D_{50,t}$ forecast was 43% as opposed to 47% for the null model, indicating that the model-averaged forecast model was not effective at reducing the average magnitude in $\widehat{\text{EOS}}_{d,t}$ prediction errors. Additionally, the model-averaged forecast was ineffective at reducing the positive bias in these predictions from the null model as indicated by the MPE for both models being approximately 10% (Table \@ref(tab:eos-table)). 

A visual example of prediction accuracy and uncertainty from two recent years is provided in Figure \@ref(fig:eos-preds). The upper panels show the time series of $\widehat{\text{EOS}}_{d}$ when the null and model-averaged forecast models were used to inform the location of the logistic cumulative timing curve in 2013. The horizontal line shows the observed value of $\text{EOS}_t$. 2013 is an example of when the null model would have been preferable to use (in terms of accuracy) and 2014 shows a case when the model-averaged forecast would have performed better. According to Figure \@ref(fig:forecasts), 2014 was one of the earliest runs on record, which explains using $D_{50}$ informed by the null model lead to over estimates of $\widehat{\text{EOS}}_{d}$ for much of the season that year.

### Run timing _versus_ run size relationship {#rt-n-results}

\noindent
There appeared to be little evidence to lend support for the hypothesis that run timing and run size are related for the Kuskokwim River Chinook salmon stock. Based on the visual depiction of the relationship (Figure \@ref(fig:rt-n)), there appeared to be a weak negative pattern. However, the fitted regression model suggested the effect of run size on run timing was weak: on average, $D_{50}$ occurred 1.15 (95% confidence limits; -2.58 -- 0.27) days earlier for each 100,000 fish increase in run size, which was not significantly different than no effect of run size on run timing ($p = 0.11, R^2 = 0.08, \hat{\sigma} = 3.6$). Additionally, based on forecast cross-validation, the model that included run size did not perform better at prediction than the null model. On average, the model that included run size model resulted in an estimated absolute forecast error of 0.3 (95% confidence limits; 0.47 -- 1.07) days larger than that of the null model ($p = 0.42$).

## Discussion

\noindent
The environmental relationships with run timing I detected for the Kuskokwim River Chinook salmon stock are consistent with patterns found elsewhere in the region [_e_._g_., @mundy-evenson-2011; @hodgson-etal-2006]. Specifically, I found that warmer years were typically associated with earlier-than-average runs as were years with lower-than-average SIC. These findings are consistent with the water column stability hypothesis suggested by @mundy-evenson-2011. The amount of unexplained variation in the Kuskokwim model appears to be comparable with the Yukon River Chinook salmon stock as well [@mundy-evenson-2011]. Using the Kuskokwim River relationships shown in Figure \@ref(fig:relationships), the correlation with $D_{50}$ was -0.52, -0.57, and 0.59 for air temperature, SST, and SIC, respectively. For the Yukon River Chinook salmon stock, @mundy-evenson-2011 found correlations of -0.59, -0.72, and 0.66 for the same variables but measured at different spatial and temporal scales and with approximately 10 more years of data included [Table 2 in @mundy-evenson-2011]. These similar correlations indicate the signals given by environmental variables are of relatively equal strength between these two systems.

Given the presence of the environmental relationships, it is somewhat surprising that the null model forecast performed better on average than did the model-averaged forecast. This could, potentially, be due to the fact that a variety of biological [size; _e_._g_., @bromaghin-2005 and morphology; @hamon-etal-2000] and abiotic factors [temperature; _e_._g_., @salinger-anderson-2006, river discharge; _e_._g_., @keefer-etal-2004, and migration distance; _e_._g_., @eiler-etal-2015] may affect migration rate (and subsequently, encounter probability) and catchability, introducing additional variability in my run timing estimates. Future research that accounts for these effects on encounter probability or catchability could offer improved predictions of run timing. Regardless of the underlying drivers, the overall prevalence of years with average run timing likely led to the enhanced performance of the null model.

Although the null model performed better in the long-term average (_i_._e_., lower $\overline{\text{AE}}$ as of 2018), there are reasons a manager may still justifiably prefer the model-averaged forecast. First, the difference in $\overline{\text{AE}}$ between the model-averaged forecast and the null model was 0.5 days, which is small relative to the amount of annual variation in run timing (a 17 day range for $D_{50}$ over 35 years). Second, the model-averaged forecast performed equally well in terms of forecast accuracy regardless of the type of run timing it was used to forecast (_i_._e_., prediction error equal in extreme early/late and average years; Figure \@ref(fig:mae-subsets)b). In contrast, the null model only performed comparably well in years with run timing within $\pm 3$ days from average and error increased precipitously in more extreme years. Third, the 95% prediction intervals from the null model seemed too wide as 100% of the observations fell within the intervals, whereas 92% of the observations fell within the prediction intervals from the model-averaged forecast (which is closer to the ideal coverage, _i_._e_., 95%). Prediction uncertainty was lower under the model-averaged forecast than the null model, which could ultimately lead to fewer run timing scenarios being considered early in the season.

As for the value of having access to a run timing forecast to in-season run size assessments, I showed that the model-averaged forecast offered no greater performance in terms of accuracy than the null model. However, the key difference between approaches was the reduced uncertainty in $\text{EOS}_t$ predictions when using the model-averaged forecast due to the exclusion of extreme early or late runs which lead to extreme low and high $\text{EOS}_t$ predictions early in the season. The null model was forced to always consider these scenarios, resulting in greater uncertainty (coefficient of variation) in $\text{EOS}_t$ predictions, particularly between June 15 and June 30 when key management decisions are made. Due to the large amount of uncertainty under the null model (which is essentially the currently used method), $\text{EOS}_t$ predictions go largely ignored for much of the season and the pre-season run size forecast is trusted instead. If the environmental variable forecast model were to be used, it likely would provide managers with more information when making decisions. It should be emphasized that I did not evaluate the ability to assess actual run size herein, only an index of run size ($\text{EOS}_t$). Though because $\text{EOS}_t$ can be expressed as a function of actual run size (using assumptions about catchability), more precise predictions of $\text{EOS}_t$ should presumably result in more precise predictions of actual run size.

The sliding window algorithm was an objective, adaptive, and data-driven tool for temporal selection of environmental predictors and therefore may be more appropriate than choosing a time window based on _a priori_ assumptions, particularly in forecasting applications. The algorithm relied on the data and predictive performance to select the window used for the next year's forecast. This framework allowed for predictor variables to change adaptively as a more accurate window became apparent. This quality of the sliding window algorithm makes it intuitive and potentially preferable in the face of a changing climate. However, the algorithm was computationally intensive for the retrospective analysis. The R code to select windows for all variables/grids for the procedure took approximately 1.5 days to complete on a desktop computer with a 3.60 GHz processor with four cores and 32 GB of RAM. Each year took approximately 3.5 hours to complete (depending on the number of years the forecast cross-validation procedure was conducted on). The flexibility of the approach also hinders the ability to typify a year as "warm" or "cold", as these criteria may change when an additional year of data is included and a new window is selected. An additional drawback of the sliding window algorithm is that it may be difficult to explain to managers and stakeholders, which may lead to confusion and distrust in the method.

Model-averaging across grid cells for SST and SIC was also an objective, adaptive, and data-driven, (though computationally intensive) solution for dealing with the spatial nature of predicting run timing from these two variables. Of course, it would be possible to average all daily values across grid cells each year and perform the sliding window algorithm on these means. However, this would ignore the fact that some grid cells inherently have a stronger timing signal and would likely insert more variation into the predictive relationship. Additionally, the model had the flexibility to place more weight on different grid cells as more years of data were added, again adding to the flexibility of our overall approach which may be preferable in the face of a changing climate. However, the inherent complexity of including the spatial structure again makes it more difficult to typify a year as "warm" or "cold" as there are many values each year for SIC and SST and the strength of the run timing signal given by each grid cell varies.

These two complexities to my analysis (sliding window selection and model-averaging across spatial grid cells) made the interpretation of effect sizes and selected windows difficult in a biologically meaningful way because the best windows and spatial grid cells could change from year to year. I generally see this flexibility based on predictive performance as more important for this particular analysis than biological inference on research topics like determining the most influential variable on run timing variation or determining the migration route of Chinook salmon through Kuskokwim Bay. These examples remain exciting research questions for the future, however I see them as secondary objectives of importance. My focus was on the prediction of a single critical quantity, $D_{50}$, which could aid in-season decision-making. A separate issue confounding biological interpretation is that each variable had some direct or indirect link to temperature, suggesting that there is strong potential for multicollinearity among predictor variables. It is well known that correlated predictor variables can result in biased coefficient estimates and variance inflation [@neter-etal-1996]. This was one reason coefficient estimates were only presented in the single-predictor case (Table \@ref(tab:coefs-table)), as I caution against their interpretation in this particular case. However, my focus was entirely on predictive ability, which is generally known to be unaffected by multicollinearity [@graham-2003].

An important caveat of this analysis is that I assumed a negligible influence of downstream harvest on the ability of the BTF to index the true Chinook salmon run timing. This assumption was likely violated in some years but the magnitude of the impact is unknown. From 2008 to 2017, the approximate average exploitation rate by only villages downstream of the BTF index was 20% (_versus_ 38% for villages across the whole drainage), thus there is the potential for a downstream harvest bias on perceived run timing (_i_._e_., approximately a fifth of the run can be removed before it is sampled by the BTF). Moderate to high exploitation rates would not necessarily bias the BTF index if the timing of the harvest was similar to the run. However, in the Kuskokwim, subsistence harvest has historically focused on the early portion of the run [@hamazaki-2008]. When coupled with moderate exploitation rates (_i_._e_., 20%), this early nature of the fishery is likely to have resulted in detected timing curves that were biased late (due to early fish being removed before they are sampled by the BTF) to some unknown degree. Historical and future interpretation of the BTF is further complicated by the operation of the fishery, such as a recent regulatory measure which mandates that no directed Chinook fishery may begin on or before 11 June. I suspect that the magnitude of the bias in the index due to the timing of downstream harvest would be small and would not likely affect the general conclusions of this analysis (although residual variation in environmental-run timing relationships would likely be lower if accounted for, which could improve forecast performance). I suggest that future studies should attempt to develop methods that remove harvest effects from the BTF index and other similar indices and assess the magnitude of potential bias. Even if harvest bias could be removed from the historical index values, addressing bias in the test fishery index would be unfeasible during the season because spatially and temporally explicit harvest data are often unavailable until the season has concluded, and the data regarding the temporal distribution are fragmentary.

It was unsurprising that no meaningful relationship exists between estimated run size and run timing. Given that small/early and large/late runs are problematic for in-season management [@adkison-cunningham-2015], I see the lack of a relationship as beneficial to the management effort. In other words, a small run is no more likely to be early than it is to be late, and the same is true of large runs. For managers, this means that although these small/early or large/late scenarios have occurred in the past, they need not be particularly worried about them due to an overwhelming prevalence over other run size/run timing scenarios.

There is evidence to suggest that, on a population demographic scale, substock structure and relative stock composition may influence the run timing of the aggregate. For example, @clark-etal-2015 showed that Chinook salmon that travel farther in the drainage to spawn (_i_._e_., headwaters) enter the main-stem earlier in the season. This point is supported in the Kuskokwim River based on ADF&G radio telemetry data [@smith-liller-2017a; @smith-liller-2017b; @stuby-2007], which showed that the date at which 50% of headwaters fish were tagged occurred as many as 10 or 11 days earlier than tagging of fish bound for middle river and lower river tributaries. Thus, it is more appropriate to view the timing curve detected by the BTF index as a mixture distribution made up of several distinct substocks, each entering at different times. The cumulative effect of this is one curve that looks logistic likely because the various substocks overlap to a large extent. However, it is not difficult to see that if in some years the headwaters substocks made up a greater proportion of the aggregate stock than the lower and middle river substocks, the timing curve of the aggregate would be earlier than if other stocks had a greater contribution. Using genetic techniques, @anderson-beer-2009 found that variations in the relative abundances of the populations composing the spring Chinook salmon run in the Columbia River, USA, explained 62% of the variation in annual run timing. This is a source of variation that was not accounted for in this analysis for at least two reasons. First, the resolution to divide the aggregate curve into its substock components is not currently available: Kuskokwim River Chinook telemetry studies were conducted from 2003 -- 2007 and 2015 -- 2017 and the aggregate timing curve does not deviate enough from the smooth logistic curve to separate the different substock components. Second, information on the relative contribution both in the past and in the forecast year would be necessary to include this complexity in the run timing forecast. This detailed level of substock information is not available for the Kuskokwim River. The telemetry data can shed some light on these issues, but they are confounded by factors like harvest timing [some stocks may be harvested preferentially purely due to the timing of the fishery, which does not mirror that of the aggregate run; @hamazaki-2008] or the potential of tagging stock components in some proportion other than their true contribution.

Methods exist to incorporate run timing forecasts from analyses like these into in-season assessment and management efforts. My predictions of $\widehat{\text{EOS}}_{t,d}$, which is an index of run size, could be used to predict the total end-of-season run size on each day using a regression model that relates historical reconstructed total run abundance and observed $\text{EOS}_t$. These in-season run predictions could be used to update pre-season run size forecasts with in-season data using methods such as inverse variance weighting [_e_._g_., @walters-buckingham-1975] or Bayesian inference [_e_._g_., @fried-hilborn-1988]. Information updating may be preferable in cases when the pre-season run forecast is biased, because it would allow for the perception of run size to pull away from the forecast when in-season data suggest it is highly unlikely. As I have shown here, uncertainty in $\widehat{\text{EOS}}_{t,d}$ predictions is a function of the precision in the anticipated proportion of the run completed-to-date \@ref(eq:eos). My analysis suggests that incorporating run timing forecasts into estimates of $p_{d,t}$ (and thus $\text{EOS}_t$) may provide managers with more certainty regarding interpretation of in-season abundance indices, which would facilitate updating of pre-season forecasts with data from the run. 

Since completion of the run timing forecast model, I have had the opportunity to more thoroughly test the performance of using it to estimate run size in-season in a Bayesian framework [@staton-catalano-2019]. We found that the inclusion of information from the timing forecast model had no influence on the predicted run size, and thus that method had the same performance as a method that simply expected years with average run timing. This finding was a result of the large amount of sampling variability in the BTF: even when the $\text{EOS}_{t}$ value is completely observed in year $t$, posterior abundance predictions are still highly uncertain (coefficient of variation of 20%) as a result of a scattered historical relationship. The main conclusion here is that although the migration timing forecast had utility for reducing uncertainty of predictions for $\widehat{\text{EOS}}_t$, $\text{EOS}_t$ appears to be a generally poor predictor of run abundance. This indicates that if a more reliable index (such as a sonar project) with a tighter relationship with the total run size were to become available, the run timing forecast model developed here may have greater utility.

\clearpage

\begin{singlespace}

```{r rt-ests-table}
ests = read.csv("tab/Ch2/MLE_run_timing_estimates.csv")

d50 = ests[ests$param == "d50",]
d50[,3:4] = round(d50[,3:4], 1)
h = ests[ests$param == "h",]
h[,3:4] = round(h[,3:4], 3)
paren = function(x) {
  paste(x[1], " (", x[2], ")", sep = "")
}

tab = data.frame(
  year = sort(unique(ests$year)),
  d50 = apply(d50[,3:4], 1, function(x) paren(as.numeric(x))),
  h = apply(h[,3:4], 1, function(x) paren(as.numeric(x)))
)

colnames(tab) = c("Year", "$D_{50,t}$", "$h_t$")

kable(tab, "latex", longtable = F, booktabs = T, escape = F, align = "ccc", linesep = "",
      caption = "Parameter estimates (mean with standard error in parentheses) from
      logistic curves given by (\\ref{eq:logistic}) fitted to each year separately. $D_{50,t}$ is expressed
      as the day-of-the-year, for reference, day 174 is June 22 in a leap year and June 23 in a normal year.") %>%
  kable_styling(latex_options = c("repeat_header"))
```

\end{singlespace}

\clearpage

```{r scwa-dates-table}
tab = data.frame(
  var = c("AIR", "PDO", "SST", "SIC"),
  dfdoy = c(1,"$\\--$",92,50),
  dfdate1 = c("Jan. 1", "Jan.", "Apr. 2", "Feb. 19"),
  dfdate2 = c("Jan. 1", "$\\--$", "Apr. 1", "Feb. 19"),
  dldoy = c(151,"$\\--$",151,130),
  dldate1 = c("May 31", "Jun.", "May 31", "May 10"),
  dldate2 = c("May 30", "$\\--$", "May 30", "May 9")
)

kable(tab, "latex", booktabs = T, align = "lcccccc",
      caption = "The input constraints used in the SCWA (Section \\ref{scwa}) for selecting the best time period of prediction for each covariate. Note that only monthly variables were available for PDO.",
      col.names = linebreak(c("Variable", "DOY", "Non-Leap\nYear", "Leap\nYear", "DOY", "Non-Leap\nYear", "Leap\nYear"), align = "c"), escape = F) %>%
  add_header_above(c(" " = 1, "$\\\\boldsymbol{D_0}$" = 3, "$\\\\boldsymbol{D_n}$" = 3), escape = F) %>%
  column_spec(1, bold = T)
```

\clearpage

```{r coefs-table}
tab = read.csv("tab/Ch2/Single Variable Regression Estimates.csv", stringsAsFactors = F)
model = toupper(tab$X)
tab$t = tab$slope/tab$slope.se
tab = tab[,c("intercept", "slope", "t", "r.sqr", "sigma", "f.stat")]
tab = round(tab, 2)
# tab = cbind(tab, c("*", "*", "**", "**"))
tab = cbind(model, tab)
tab[,"f.stat"] = paste(tab[,"f.stat"], footnote_marker_symbol(c(1,1,2,2), format = "latex"), sep = "")

colnames(tab)[1] = "Variable"
colnames(tab)[2] = "$\\boldsymbol{\\hat{\\beta}_0}$"
colnames(tab)[3] = "$\\boldsymbol{\\hat{\\beta}_1}$"
colnames(tab)[4] = "$\\boldsymbol{t}$"
colnames(tab)[5] = "$\\boldsymbol{R^2}$"
colnames(tab)[6] = "$\\boldsymbol{\\hat{\\sigma}}$"
colnames(tab)[7] = "$\\boldsymbol{F}$"

kable(tab, "latex", booktabs = T, escape = F,
             caption = "Estimates and statistics of the effects of each of the four single-variable forecast models fitted with all  $D_{50}$ and environmental data through 2018.", align = c("c")) %>%
  kable_styling(position = "center") %>%
  column_spec(column = 1, bold = T) %>%
  footnote(
    symbol_title = "Significance Codes",
    symbol = c("P < 0.01", "P < 0.001"))
```

\clearpage

```{r eos-table}
tab = read.csv("tab/Ch2/eos_pred_table.csv", stringsAsFactors = F)
dates = tab$date
rnded = round(tab[,-1], 2) * 100

tab = apply(rnded, 2, function(x) paste(x, "%", sep = "\\"))
tab = cbind(dates, tab)

colnames(tab) = c("\\textbf{Date}", "\\textbf{NULL}", "\\textbf{Mod. Avg.}", "\\textbf{NULL}", "\\textbf{Mod. Avg.}", "\\textbf{NULL}", "\\textbf{Mod. Avg.}")

kable(tab, "latex", booktabs = T, align = "ccccc",
      caption = "Retrospective accuracy and uncertainty of predictions of the end-of-season cumulative CPE at the Bethel Test Fishery ($\\widehat{\\text{EOS}}_{d,t,i}$) as informed using two methods of obtaining estimates of the fraction of the run complete ($p_{d,t}$) as described in the text (Section \\ref{eos-preds}). Accuracy is expressed as the mean percent absolute percent error (MAPE) and mean percent error (MPE) and uncertainty in the prediction made each day and year is expressed using the coefficient of variation (CV).", escape = F) %>%
  add_header_above(header = c("", "MAPE" = 2, "MPE" = 2, "CV" = 2), bold = T)%>%
  column_spec(1, bold = T)
```

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/map.png}
  \caption{Map of Kuskokwim Bay where Chinook salmon likely stage for transition to freshwater. Shows grid cells from which daily SST values were used. Daily SIC values came from the same grid cells, though excluding grid cell 45 due to missing values.}
  \label{fig:ch2-map}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/p-ccpue.png}
  \caption{Shape and variability of run timing patterns of the Kuskokwim River Chinook salmon stock as sampled by the Bethel Test Fishery, 1984 -- 2018. Each grey curve represents a year standardized by the total end-of-season cumulative CPE and the black line represents the average value across years on each day of the season.}
  \label{fig:p-ccpue}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/relationships.png}
  \caption{Relationships between the four single environmental variables and run timing $\left(D_{50}\right)$ using data from optimal climate windows when 2018 was included in the training data. For illustration purposes only, gridded variables SST and SIC were combined by weighted averaging, where the weight of each grid cell was assigned the $\text{AIC}_{\text{c}}$ weight of that grid cell when grid cell-specific models were fit. Grey bands are 95$\%$ prediction intervals.}
  \label{fig:relationships}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/window-changes.png}
  \caption{Changes in selected climate windows as more years of training data were added in the retrospective forecasting analysis. Bottom and top lines show the first and last day of the selected climate window, respectively, as more years were added. The year axis corresponds to the selected window after including environmental and run timing data from that year in the training data. For example, the windows shown for 2017 were used to produce the forecast for 2018. Panel (\textit{a}) is Bethel air temperature, panels \textit{b}1 -- \textit{b}4 are SST windows for four sample grid cells and panels \textit{c}1 -- \textit{c}4 are SIC windows for the same four sample grid cells. Sample grid cells from Figure \ref{fig:ch2-map} shown for SST and SIC are as follows: grid cell 8 (\textit{b}1, \textit{c}1), grid cell 44 (\textit{b}2, \textit{c}2), grid cell 12 (\textit{b}3, \textit{c}3), and grid cell 48 (\textit{b}4, \textit{c}4). Selected windows for PDO are not shown because the single month of May was selected in all years.}
  \label{fig:window-changes}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/forecasts.png}
  \caption{Observed and forecasted median run date ($D_{50}$) under the three approaches: the null (intercept only/historical average) model, the single regression model with the lowest forecast cross-validation score up to the forecasting year, and a model-averaged forecast that combined forecasts from 16 nested regression models. Black points/lines are the time series of $D_{50}$ detected by the BTF. Grey points are out-of-sample forecasts with 95$\%$ prediction intervals shown as error bars. $\overline{\text{AE}}$ and $\widetilde{\text{AE}}$ are the mean and median absolute forecast errors from 1995 to 2018, respectively.}
  \label{fig:forecasts}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/ae-changes.png}
  \caption{Evolution of mean ($\overline{\text{AE}}$) and median ($\overline{\text{AE}}$) absolute forecast error under the three investigated forecasting approaches: the null (intercept only/historical average) model, the single regression model with the lowest forecast cross-validation score up to the forecasting year, and a model-averaged forecast that combined forecasts from 16 nested regression models. Each point is the average of absolute errors of all years before and including the corresponding year on the $x$-axis, starting in 1995, which was the first year out-of-sample forecasts were evaluated.}
  \label{fig:ae-changes}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/mae-subsets.png}
  \caption{$\overline{\text{AE}}$ under three forecast approaches calculated by either (\textit{a}) including years with a $D_{50}$ value within $\pm x$  days of the all-year average or (\textit{b}) including years with a $D_{50}$ value outside $\pm x$ days of average, where $x$ is the number of days indicated on the $x$-axis. Bottom panels show the number of observed years in which the appropriate $\pm x$ days criterion was met. Shaded regions in the hypothetical distributions show the types of $D_{50}$ values that were included in the calculation of $\overline{\text{AE}}$. One point that may enrich inference from this figure (and is shown in the shaded normal distributions) is that panel (\textit{a}) becomes more inclusive from left to right by adding years that are more dissimilar to the average in the calculation of $\overline{\text{AE}}$ whereas panel (\textit{b}) becomes more exclusive from left to right by removing years that are similar to the average.}
  \label{fig:mae-subsets}
\end{figure}

\clearpage

\begin{figure}
  \centering
  \includegraphics{img/Ch2/eos-preds.png}
  \caption{In-season predictions of end-of-season cumulative BTF CPE under the model-averaged forecast using environmental variables and the forecast under the null model in 2013 (panels \textit{a} and \textit{b}) and 2014 (panels \textit{c} and \textit{d}). Intended to illustrate cases in which a manager would benefit from having access to the model-averaged run timing forecast model using environmental variables (2014) and when the null model would have performed better (2013). Horizontal lines are the true end-of-season cumulative BTF CPE and light grey regions are 80$\%$ confidence intervals obtained \textit{via} Monte Carlo simulation. Dark grey vertical regions indicate the period when key harvest decisions are made with respect to Chinook salmon in the Kuskokwim River.}
  \label{fig:eos-preds}
\end{figure}

\begin{figure}
  \centering
  \includegraphics{img/Ch2/rt-n.png}
  \caption{Relationship between $D_{50}$ and run size for Kuskokwim River Chinook salmon with two fitted models shown: the null model (which assumed constant mean $D_{50}$) and the run size model (which assumed the mean $D_{50}$ changes as a function of run size). As described in the text (Section \ref{rt-n-results}), the effect of run size on run timing was weak and not significantly different than no effect. Additionally, knowledge of run size did not result in smaller average prediction errors of $D_{50}$ than not having this knowledge.}
  \label{fig:rt-n}
\end{figure}