Many populations of Sockeye Salmon *Oncorhynchus nerka* in the eastern North Pacific Ocean experienced significant productivity declines that began about 1990, but there is no consensus on the mechanisms responsible. To better understand Sockeye Salmon survival trends, we examined the 50-year time series for two age-classes of Sockeye Salmon smolts from Chilko Lake in central British Columbia. Arranging survival time series for both age-classes by ocean entry year and combining them, weighted by a proxy model of sampling variance, reduced the sampling variance in the original age-1 smolt survivals sufficiently to indicate a linear trend of increasing survival from 1960 to 1990 that suddenly changed at or near 1991 to a lower and declining trend from 1992 to 2008. Neither density nor mean length influenced smolt survival. Returns in a given year were not good predictors of siblings returning in subsequent years. Time spent at sea increased linearly beginning around 1970. Although smolt survivals differed between ecosystem regimes, there was only the one clear pattern break about 1991. To improve our understanding of mechanisms, survival trends were compared with environmental indices that included catches and hatchery releases of potentially competing salmon from around the North Pacific Ocean. Smolt survivals were more similar to abundance indices of Sockeye Salmon, Chum Salmon *O. keta*, and Pink Salmon *O. gorbuscha* than to indices of global, regional, or local ocean climate. Our results are consistent with the hypothesis that salmon productivity in the North Pacific declined soon after 1990. We present a simple model to illustrate how increased competition at sea, related to the release of large numbers of hatchery salmon, in conjunction with changes in ocean productivity, may have played a significant role in improving Sockeye Salmon survivals while reducing their growth before 1991. After 1991, these factors may have acted to reduce survivals while the growth of survivors showed no effect.

Anadromous Pacific salmon *Oncorhynchus* spp. comprise a multispecies complex of varying productivities. Their recent abundance in the Pacific Ocean, as reflected by commercial catch, is as high as it has ever been (Irvine and Fukuwaka 2011). Global abundances are driven primarily by Pink Salmon *O. gorbuscha* and Chum Salmon *O. keta*, as well as, particularly in the eastern North Pacific Ocean, by Sockeye Salmon *O. nerka* (Eggers 2009; Ruggerone et al. 2010; Irvine and Fukuwaka 2011). The status of Sockeye Salmon populations varies among regions however, and in British Columbia's Fraser River, low numbers of returning salmon in recent years are a major concern (Grant et al. 2011; Rand et al. 2012).

The Fraser River watershed is one of the world's greatest salmon producers (Northcote and Larkin 1989), although numbers returning annually are highly variable. Sockeye Salmon are the most economically valuable salmon species in the watershed, and have provided a commercial harvest since the early 1870s (Meggs 1991) and a First Nations (native North American) harvest for millennia. Returning Sockeye Salmon are divided into three major groups, based on run timing, that comprise at least 22 Conservation Units, biological groups of salmon that are genetically or ecologically distinct from each other (Holtby and Ciruna 2007). Declining Sockeye Salmon returns from 1992 to 2009, and an exceptionally low 2009 return (smallest since 1947) resulted in the Canadian government establishing a judicial inquiry on Fraser River Sockeye Salmon (Cohen 2012a, 2012b, 2012c). Ironically, the inquiry was barely underway when 2010 saw the largest Sockeye Salmon return in the previous 100 years. The inquiry represented 2.5 years of work that included 128 d of evidentiary hearings, 2,145 exhibits, and testimony from 179 expert witnesses (Cohen 2012c:86). Fifteen detailed technical reports plus various primary publications (e.g., Beamish et al. 2012; Connors et al. 2012; Peterman and Dorner 2012; Preikshot et al. 2012; Thomson et al. 2012) were produced. Beamish et al. (2012) attributed the low returns for Fraser River Sockeye Salmon in 2009 to local effects within the Strait of Georgia, while Thomson et al. (2012) documented the importance of various factors at different life history stages in determining survival. The inquiry concluded that marine factors in particular were implicated in the broad-based regional decline of Sockeye Salmon stocks, but was unable to determine the relative importance of specific factors (Cohen 2012c:88).

Temporal patterns in salmon survival have been influenced by many factors including major ecosystem regime shifts in 1977 and 1989 (e.g., Beamish and Bouillon 1993; Hare and Mantua 2000; Irvine and Fukuwaka 2011). Some researchers also identified 1999 as a shift (e.g., Peterson and Schwing 2003). With respect to Fraser River Sockeye Salmon, Beamish et al. (2004b) found survivals were high following the 1977 shift and declined after 1989. In contrast, McKinnell and Reichardt (2012) reported no increase in Fraser River Sockeye Salmon survival after 1977, although they found first-year marine growth abruptly declined about 1977. Ruggerone et al. (2003) attributed reduced survival of Sockeye Salmon in Alaska to increased competition with Pink Salmon. In the case of Fraser River Sockeye Salmon, Connors et al. (2012) concluded the effect of competition with Pink Salmon could be exacerbated if Sockeye Salmon were exposed to farmed salmon during their out-migration, particularly during warm ocean conditions. Ruggerone et al. (2010) and others have raised concerns that density-dependent interactions in the ocean resulting from hatchery releases may reduce Pacific salmon growth and survival.

Many analyses of Pacific salmon survivals rely on stockrecruit data. Results from these studies, often conducted across many stocks, have greatly improved our understanding of the underlying patterns of Sockeye Salmon growth (Peterman 1984) and survival (Peterman et al. 1998; Pyper et al. 2005; Peterman and Dorner 2012). For example, Peterman and Dorner (2012) demonstrated productivity declines since the 1990s or earlier for many Sockeye Salmon stocks from Washington to southeast Alaska and suggested that climate-driven increases in mortality induced by pathogens, as well as increased predation or reduced food due to oceanographic changes, may be the explanation. Yet stock-recruit data typically do not allow one to separate effects occurring at different life history stages. Of the numerous Fraser River Sockeye Salmon populations, there are only two, Cultus Lake in the lower Fraser River watershed and Chilko Lake in the central, where survival is estimated before and after smolts leave their rearing lake. The Cultus Lake time series is relatively short, but at Chilko Lake, 775 km upstream from the ocean, a reasonably consistent and nearly unbroken time series of abundance and length by age for spawners and smolts has been collected since the early 1950s. Most Sockeye Salmon smolts exit Chilko Lake during spring at age 1, but some spend a second year in freshwater. Age-1 and age-2 smolts leaving the lake during the same year have different parents and a different freshwater life history but are exposed to similar marine conditions. Previous researchers examining the Chilko Sockeye Salmon time series generally focused on the predominant life history, age-1 smolts, most of which return as adults after 2.5 years at sea (e.g., Henderson and Cass 1991; McKinnell 2008; Grant et al. 2011). Although survival estimates for age-2 smolts have a larger sampling error than estimates for the more abundant age-1 smolts (Bradford et al. 2000), we recognized that information gleaned from age-2 smolt survivals could help to clarify survival patterns and understand the processes shaping survival and population dynamics for Sockeye Salmon.

To better understand the role of marine factors affecting Sockeye Salmon survival, we focused on the postsmolt life history of Chilko Lake Sockeye Salmon. We compared age-1 smolt survival estimates with those of age-2 smolts leaving the lake in the same year, investigated uncertainty associated with these estimates, and developed a new survival time series. We evaluated the influence on survival of smolt abundance (density dependence), smolt quality (age, length, and condition), and environmental trends (various North Pacific and regional ocean climate indices). Finally, we examined the potential effects on Sockeye Salmon survival and growth patterns of density-dependent processes at sea as represented by indices of multispecies salmon abundance.

## METHODS

*Sockeye Salmon life history and data sources.*—A brief synopsis of the life history of Chilko Lake Sockeye Salmon and how data were gathered follows. We began our analyses with the fish that were spawned in 1958, the beginning of a continuous series. Smolt numbers were not estimated in 1991. Smolt and spawner sampling methodologies were reviewed by Henderson and Cass (1991), Roos (1991), and Grant et al. (2011).

Chilko Lake Sockeye Salmon start life as fertilized eggs in the fall of the brood year, emerge as fry from spawning sites the next spring, and spend 1 or 2 years in Chilko Lake before exiting as smolts during spring. When smolts exit the lake after usually 1 year, they pass through an enumeration fence at the lake outlet where their numbers are estimated photographically. A random sample of smolts is measured daily for length, and smolts larger than a specified length are collected for age determination. That specified length, typically but not always 85 mm (K. Brenner, Fisheries and Oceans Canada [DFO], Kamloops, personal communication) is smaller than the minimum length of age-2 smolts. The first smolt records we analyzed were collected in 1960, recording age-1 smolts (*S*_{1}) from brood (spawning) year 1958 and age-2 smolts (*S*_{2}) from brood year 1957. In addition, samples of age-1 and age-2 smolts (determined by size) were taken at the counting facility, from near the beginning, the peak, and the end of most smolt runs, often several hundred smolts per year, and preserved in 10% formalin. Subsequently, ages were confirmed from scales, and fish were weighed and measured for length to evaluate changes in fish size and condition (J. Hume, DFO Cultus Lake, personal communication).

Most *S*_{1} return to freshwater as maturing adults in late summer after two winters in the ocean and were categorized as *R*_{1.2,} meaning they lived a total of three winters, one in freshwater (in addition to one winter as a developing egg), and two in the ocean (Figure 1). A small proportion of *S*_{1} spend one or three winters at sea, returning as *R*_{1.1} and *R*_{1.3}. Age-2 smolts (*S*_{2}) stay for a second growing season and a second winter in Chilko Lake and also typically spend two winters in the ocean (*R*
_{2.2}), but again some return after one or three winters at sea, as *R*_{2.1} and * R_{2.3}.* Some life history stages such as

*S*

_{3}and

*R*

_{2.3}were rare and we have no records of

*R*_{3.2}.The Pacific salmon fisheries were sampled for stock and age composition to estimate the numbers of Chilko Lake Sockeye Salmon mixed with other salmon stocks in catches, on their return migration. Allocation of Sockeye Salmon to discrete stocks was based on distinctive freshwater growth patterns seen on the scales (Gable and Cox-Rogers 1993), supplemented by DNA markers beginning in 2000 (Grant et al. 2011).

The number of fish escaping fisheries and returning to Chilko Lake to spawn were estimated prior to 2009 by mark-recapture (e.g., Schubert and Fanos 1997), except for 1967 when they were estimated by expanding visual counts at Henry's Bridge, 12 km below the lake (Grant et al. 2011). Sockeye Salmon carcasses were sampled throughout the Chilko Lake spawning areas to obtain scales and measure length. Scales and otoliths were collected so that estimates of spawner abundance could be partitioned into age-classes. Since 2006, imaging sonar (DIDSON) has been used to estimate spawner abundance (Cronkite et al. 2006; Holmes et al. 2006). Numbers estimated by DIDSON and by mark-recapture were similar and, since 2009, DIDSON estimates have replaced mark-recapture (Benner, and T. Cone, DFO Annacis Island, personal communication[s]). Total adult returns, as catch in the fisheries plus escapement to the spawning grounds, were estimated each year and then, based on the ages of spawners collected at Chilko Lake, divided into six life history categories (Figure 1).

Terms and abbreviations for this analysis are described in Table 1. Variables are in italic uppercase: e.g., *S*_{
1,t
}
*, R_{1.2,t+2}, SAR_{1}.* Fitted regression parameters are in italic lowercase: e.g.,

*a*,

*b*Means were reported with the standard error in parentheses, i.e., (SE), unless explicitly stated otherwise. Correlation probabilities were corrected for shared autocorrelation by estimating the effective degrees of freedom (Pyper and Peterman 1998).

_{1}.
*Smolt survival estimates.*—Smolt survival was referenced to the ocean entry year t when smolts left Chilko Lake and were first observed. This survival calculation included the period when smolts migrated downstream from Chilko Lake to the ocean, as well as their time at sea. The terms *S*_{1,t}, and *S*_{2,t+ 1} thus refer to age-1 and age-2 smolts from the same brood year but having different ocean entry years. The total returns from age-1 smolts,

were the survivors of *S*_{1,t} from brood year t — 2 that entered the ocean in year *t* and returned 1, 2, or 3 years later. Similarly, * R_{2,t}
* were the survivors of

*S*

_{2,+ 1}from brood year t — 3 that also entered the ocean in year

*t.*The estimate of postlake survival for

*S*

_{1}is the ratio of smolts to adult returns, i.e.,

*SAR*

_{1}=*R*_{1}/*S*_{1}; SAR_{1}was a ratio of estimated variables, each with considerable sampling variance, and had a skewed distribution.

*Losses between returns and spawning.*—A Chilko Lake Sockeye Salmon that recruited to the fishery in a given year could be from one of six age-groups and four brood years (Figure 1). Recruitment for each year was rebuilt from data on abundance by age-class by brood, such that the youngest recruits in a given year are *R*_{1.1} spawned 3 years previously, and the oldest are * R_{2.3}
* spawned 6 years previously.

## TABLE 1.

Chilko Sockeye Salmon life history stages, definitions of abbreviations, and statistics for brood years 1956–2006. Spawners, smolts, recruits, and returns are expressed as millions of fish.

The time series of effective female spawners, *EFS*, was used to estimate smolt survival after fish could recruit to fisheries. Our estimate of losses between recruitment and spawning, primarily due to fisheries exploitation but including en route mortality and prespawning mortality, was

and assumed that total spawners were approximately twice the number of female spawners.

*Sampling variance.*—About 2,400 smolts and 660 spawners were sampled for age determination each year, and *S*_{2} and * R_{2}
* comprised about 4% of these samples, so the age distribution of smolts and returns had a large binomial sampling variance. There were additional sources of error in the estimates of total abundance of smolts, returns, and spawners before applying age distributions. To understand the precision of

*and*

*S*_{2}*R*

_{2}we examined the sample sizes for age data, although these were only available for recent years: 2005–2010 for smolts and 2000–2009 for returns. Binomial confidence limits (BCL) (Zar 1999:528) were calculated for smolts and returns using the R function qbinom (Pr, n,p), where n was the total number of smolts of age 1 and

*2,p =*, and Pr was 0.975 or 0.025 for the upper or lower 95% BCL. This function calculated the smallest value of

*S*_{2}/n*x*such that ƒ (x) ≥ p, where ƒ was the cumulative binomial frequency distribution (R Development Core Team 2012). Because fish below a specific size were assumed to be

*S*

_{1}and not sampled for scales, the proportion of

*S*

_{2}in the age samples were higher than in the smolt population, but we expected that BCLs for

*S*

_{2}and

*in the age samples would be proportional to those of the population and would estimate the relative precision between years for*

*R*_{2}*S*

_{2}and

*abundance.*

*R*_{2}
*Comparing and combining the* SAR_{1}*and* SAR_{2}*time* series.—We compared the time series for *SAR* _{1} and *SA R_{2}* with plots, correlations, and harmonic means (log transformations produced nearly normal distributions for survival). Principal components analysis was used to determine a common factor (PC

_{1}) from the two time series. We used generalized additive modeling (GAM, R package mgcv) on log-transformed survival as an objective comparison of the trends in

*SAR*

_{1}and

*SA*(Wood 2006; Zuur et al. 2009; R Development Core Team 2012). We tested whether a GAM smoother for each age was better than a single smoother for combined ages using the Akaike information criterion (AIC) and ANOVA.

*R*_{2}We wished to improve the time series *SAR* _{1} with information provided by *SA R_{2}
* but this required appropriate weights (w) for

*SA*because the estimates for

*R*_{2}*S*

_{2}

*and*abundances had wide BCLs in some years, which resulted in some unreliable estimates for

*R*_{2}*SA*Despite having details for age samples for onlyh a few years, we developed a model (see Appendix) to parameterize the sampling variance of

*R*_{2}.*SA*and provide appropriate weights, based on consistency between the age distribution of emigrating smolts in a year and the age distribution of returning survivors from those smolts 2 years later. These weights were used whenever the variables

*R*_{2}*S*

_{2}and

*appeared in regressions. We combined the survival estimates for each year as a weighted geometric mean that would have the same mean as*

*R*_{2}*SAR*

_{1},

where *w* was a weight from 0 to 1 and *D* was mean[log_{e}(SA*R*_{2})] — mean[log_{e}(SA*R*_{1})].

*Effect of smolt density on survival.*—We used the Ricker model, *R =* α*Se*^{βs}, to examine decreasing smolt survival with increasing smolt abundance. After log transformation,

α was estimated as exp(b_{o})- Tanaka (1962) and Larkin (1973) pointed out that independent random numbers for *R* and S in this regression will produce a significant value of *b _{1}.* McKinnell (2008) compared fits for Fraser River Sockeye Salmon stocks to equation (4) and showed that as certainty in the stock-recruit relationship decreased, estimates for α and ß increased. Incorrectly high values for Ricker model parameters suggest higher exploitation rates and lower target stock sizes, inappropriate for management under uncertainty, so a better understanding of equation (4) was important. We added the line

*f(S)*= log

_{e}[mean(R)/S] to plots from equation (4) to indicate the expected result when there was no relationship between

*R*and S, but note that

*R*approaches 0 as S approaches 0. Using nonlinear regression to fit the Ricker model indicated a bias from equation (4), but when the data were log transformed to account for increasing residuals with increasing smolt abundance, then nonlinear regression produced the same parameter estimates as did equation (4).

The second measure of smolt survival, *SA R_{2},* allowed us to develop a novel test for density dependence that avoided some of the problems of equation (4), i.e.,

where the three variables were registered to the same ocean entry year. This approach assumed that *S*_{2} and *S*_{1} coexist in the ocean such that density was effectively described by *S*_{1} (which constituted, on average, 96% of the returns of Chilko Lake Sockeye Salmon).

*Effects from ecosystem regime shifts and smolt quality.*—Our starting point for examining temporal trends in smolt survival was

where *S* was the number of smolts of age 1 and age 2 (separately) in a specific ocean entry year, *R* was the number of returns over all ages from those smolts see equation (1), and the error term, N(0,wσ), described normal residuals weighted by *w* to correct for heterogeneous variance in residuals due to heterogeneous sampling variance in * S_{2}
* and

*The x-intercept,*

*R*_{2}.*b*was the mean log survival. Because the regression slope,

_{o},*b*

_{1}

*,*was related to the correlation coefficient,

*b*

_{1}

*= r (σ*sampling variance would have resulted in

_{y}/σ_{x}),*b*

_{1}< 1 when the true value was 1.

*If b*

_{1}was significantly less than 1, indicating diminishing returns in

*R*with increasing S, we interpreted that as density-dependent survival. We preferred equation (6) to equation (4) because it avoided the ratio estimators

*SAR*

_{1}and

*SA*identified the effect of smolt abundance explicitly, and the increased range of smolt abundance reduced the effect of error in the independent variable. The residuals from this regression arose from sampling variance, environmental effects, and smolt quality, all of which were important to know.

*R*_{2},To evaluate ecosystem regime shifts, we assumed that shifts began in specific years 1977, 1989, and 1999 (as referenced in Introduction). We examined the effect on smolt survival of four regimes between 1960 and 2008 using the regression

where *P _{i}
* was a factor identifying

*i*= 1–4 regimes. We compared the fit for equation (7) to the fit for equation (6) using ANOVA to determine overall significance of the regimes, and then compared regimes using Tukey's honestly significant difference test (Zar 1999) to adjust probabilities for multiple comparisons. As a contrast to discrete regimes, we also fit a smoothly varying environmental trend to the survival from both smolt age-groups, as follows:

where *s*(*t*) was a GAM smoother for the ocean entry years. We simplified the environmental trends from equation (8) using segmented linear regression,

where *T _{i}(t)* was a factor identifying the years in different survival regimes. The result was a series of lines with intercepts

*b*and slopes

_{o,i}*b*that described temporal trends in smolt survival.

_{2,i}We added the mean annual values for smolt length, L, condition, K, and a factor for age, A_{j} to equation (9) to establish the following:

where *K* was the deviation (percent) of observed from predicted weight based on a log-log regression of preserved weights on preserved lengths (as opposed to assuming isometric growth). This was processed as forward and backward stepwise regression to assess the influence of these variables on smolt survivals.

*Variation in age composition and marine residence time.—* We examined interannual variation in the proportions of *S* _{1} and *S*_{2} that spent 1 or 3 years in the ocean, rather than the usual 2 years (Figure 1) by plotting time series of age composition of returns. We also examined correlations between returns within a brood but after differing lengths of time at sea. We removed the smolt effect from these correlations, hoping to clarify an effect due to the ocean entry year, as follows. Regression through the origin, *R*_{1.3} = bR _{1.2}, established the base case. We removed the effect of smolt numbers from *R*
_{1.3} by fitting * R_{1.3}
* =

*b*and extracting residuals which we notated as

*S*_{1}*R*

^{*}1.3. We calculated

*similarly. The regression*

*R*^{*}_{1.2}

*R*^{*}_{1.3}

**=**

*a + b*

*R*^{*}_{1.2}then tested for an effect from their common ocean entry year and from two sea winters in common.

*Marine indicators of high seas survival processes.*—We created time series plots for 24 environmental variables arranged by salmon ocean entry year and calculated the trends in these before and after 1991 for comparison to trends in *SAR.* We also examined scatter plots of *SAR* or 1*og _{e}(SAR)* against each variable and used linear regression to measure how well each variable predicted

*SAR*before and after 1991. In almost all cases, data sets covered the full duration of the

*SAR*time series. We computed results with and without the 2007

*SAR*

_{1}ocean entry year, an outlier associated with the 2009 fisheries failure, but differences were minor and we only provided results for the complete data set.

To represent density-dependent processes at sea, we used commercial catch data and hatchery release estimates from Irvine et al. (2012). Commercial catch data adequately track total salmon abundance at the species level and at the scale of the northeast and northwest Pacific Ocean (Irvine and Fukuwaka 2011). Although others (Eggers 2009; Ruggerone et al. 2010) have combined catch and spawner escapement data to develop run size biomass estimates for salmon, statistical comparisons between the time series generated by these approaches have not been made. We worked with data for Sockeye Salmon as well as for Pink and Chum Salmon since these were the most abundant salmon species and density-dependent interactions among them may be occurring (Ruggerone et al. 2010). Data were aggregated for the northeast Pacific (North America) and northwest Pacific (Asia). Catches were shifted backward by 1 year for Pink Salmon and 2 years for Chum and Sockeye Salmon to correspond to the same ocean entry year as the *SAR* index. Because interstock salmon abundance was hypothesized to be a densitydependent effect, corresponding to equation (4), we measured the ability of international salmon abundance indices to predict *log _{e}(SAR)* before and after 1991.

To evaluate large-scale ocean climate effects we used the Multivariate El Niño Southern Oscillation (ENSO) Index in May-September, the Pacific Decadal Oscillation (PDO) in November-March and May—September ( http://jisao.washington.edu/pdo/PDO.latest; Mantua and Hare 2002), the Aleutian Low Pressure Index (ALPI; Trenberth and Hurrell 1995), and the Northern Oscillation Index (NOI; Schwing et al. 2002).

To evaluate more local effects we examined the sea surface temperature (SST) data from Chrome Island in the Strait of Georgia ( www-sci.pac.dfo-mpo.gc.ca/osap/data/SearchTools/Searchlighthouse_e.htm) and the upwelling index at 48°N (ftp://orpheus.pfeg.noaa.gov/outgoing/upwell/monthly/upindex.mon). We were limited by the availability of local biological indices with time series that included the 1980s and relied upon estimates of Pacific Herring *Clupea pallasii* abundance in the Strait of Georgia (J. Schweigert, DFO Nanaimo, personal communication). Beamish et al. (2012) demonstrated that recent herring recruitments covaried with Sockeye Salmon survivals.

*Autocorrelation and heterogeneous variance.*—After examining partial autocorrelation plots for all the time series, regressions from equations (7–9) were recomputed using generalized least squares, specifically the Cor AR 1 function within the regression routine gls provided by the package nlme in the R statistical language (Zuur et al. 2009; Pinheiro et al. 2013), to estimate first-order autocorrelation and to correct regression standard errors and probabilities. Heterogeneity of variances within R2, S2, and SAR2 was addressed using the weights described in the Appendix. These weights were also used as a variance covariate via the VarConstPower function with gls. In theory, a GAM smoother, autocorrelation coefficients, and weights for heterogeneous variance can all be calculated simultaneously (e.g., GAMM, R package mgcv, R Development Core Team 2012). Zuur et al. (2009) recommended fitting a variance model first, then a predictors model. As a methodological note, we found that fitting a GAM smoother in the predictors model while simultaneously correcting for autocorrelation in the variance model resulted in these components of the regression competing for temporal patterns in variance. It was not clear how that trade-off should be resolved.

## RESULTS

### Summary Statistics and Time Series

Time series of smolts (*S*), spawners *(EFS),* and returns *(R)* varied by a factor of about 50, with numbers of *S*_{1} always exceeding *S*_{2} (Figure 2A, B) and numbers of *R*
_{1} almost always exceeding * R_{2}
* (Figure 2C, D). On average, 4.3% (SD = 4.7) of each brood of Chilko Sockeye Salmon left after a second year in freshwater (Figure 1), but occasionally the estimates for that fraction were much higher: 25.8% for brood 1965, and 16.7% for brood 1970. Age-3 smolts were only recorded recently, two each in 2009 and 2010, and were not considered in this analysis. Fish of both smolt ages usually spent 2 years at sea, but a few returned after 1 or 3 years. On average, 94.6% of the Chilko Lake Sockeye Salmon returning as adults entered the sea as

*S*

_{1}. Within that group, 1.1% were

*R*

_{1.1}(jacks) that returned to freshwater after 1 year in the ocean, 93.4% were

*that returned after 2 years, and 5.4% were*

*R*_{1.2}*R*

_{1.3}(Figure 1). Similarly for

*S*

_{2}, the proportions were 4.2%

*93.3%*

*R*_{2.1},*and 2.4%*

*R*_{2.2},*Smolt mean lengths varied by years (Table 1) but*

*R*_{2.3}.*S*

_{2}were, on average, 40% longer than

*S*

_{1}and 300% heavier. Smolt condition statistics were almost identical for both age-groups (we did not assume isometric growth).

### Precision of Age-2 Smolt Estimates

The proportion of *S*_{2} in six recent age samples ranged from 0.6% to 48% and the corresponding BCLs for *S*_{2} abundance within age samples reflected that variability. Because of lengthstratified sampling for ages, the proportion of *S*_{2} in age samples was 2–11 times greater than the proportion in the population of smolts (Table 2). For instance, the Pacific Salmon Commission (PSC) estimated that 46,940 * S_{2}
* from brood 2004 left Chilko Lake in 2007, mixed in with 77,130,000

*S*

_{1}from brood 2005. In 2007, a sample of 2,067 smolts, of which 12 were

*S*

_{2}(Table 2), was aged. The 95% BCLs for those 12

*were 6 and 19 or about ± 50% of the count in the age sample. Based on the ratio of 12*

*S*_{2}*S*

_{2}observed to 46,940

*S*

_{2}estimated, the lower and upper 95% BCLs for

*leaving Chilko Lake in 2007 were 23,471 and 74,325, respectively. This sampling variance for*

*S*_{2}*S*

_{2}amplified the sampling variance for

*SAR*

_{ 2 }

*=*because

*R*_{2}/*S*_{2}*S*

_{2}was the denominator. In 2010, as a comparison, 1,072

*S*

_{2}were observed in a sample of 3,827, or about 28%, and the 95% BCLs are only ±55, or about ±5%. The estimate for

*S*

_{2}abundance for 2010 was much more precise than the one for 2007. The PSC estimates of abundance at age involved the application of length-at-age matrices (an age—length key) to length frequency distributions, so BCLs as calculated for Table 3 only indicated the relative precision between years of abundance-at-age estimates.

## TABLE 2.

Age-2 Sockeye Salmon smolt composition of age sample and of all smolts, 2005–2010. The 95% binomial confidence limits (BCLs) for *S*_{2} in age sample are based on proportions within each age sample. The difference in the proportion of *S*_{2} in age samples compared with the proportion of *S*_{2} in estimated smolts reflects stratified sampling by length because *S*_{2} are, on average, longer than *S*_{1}.

### Precision of Adult Returns from Age-2 Smolts

For 10 recent years of available data, between 0.18% and 0.48% of the Chilko Lake Sockeye Salmon returning each year were sampled for age determination. The median sample size was about 700 (Table 3), one-third of the median sample size for smolt ages. If 20 * R_{2}
* are observed in a sample of 700, then the 95% BCLs for the number of

*that might be observed in repeated similar samples are 14 and 30. The upper confidence limit for an estimate of*

*R*_{2}*SA*would therefore be 50% higher than the estimate before considering uncertainty from

*R*_{2}*S*

_{2}and any other sources of uncertainty. If as few as seven

*R*

_{ 2 }were observed with 95% BCLs of 4 and 15, the upper confidence limit for

*SA*would be over 200% greater than the estimate. Four out of the 10 age samples had counts for

*R*_{2}*R*of less than 20 (Table 3), which suggests that roughly half of the

_{2.2}*estimates have high sampling variance (and there are additional sources of sampling variance). The sampling variance was also frequently large for*

*R*_{2.2}*and*

*R*_{1.3}*in these age samples, and the counts for*

*R*_{2.2}*were low enough (0–3) that the abundance estimates for that ageclass were overwhelmed by sampling variance. The binomial sampling variance for*

*R*_{2.3}*R*

_{1.2}was much lower than for

*The*

*R*_{2.2}.*R*

_{1.2}counted in 2009 had 95% BCLs of less than ± 3%, and the worst case in the 10 years we analyzed was 2002 at ± 5.6%.

Of 393 returns aged in 2009, three were * R_{2.2}
* from brood year 2004 and ocean entry year 2007, with 95% BCLs of 0–7. That suggested a range of ± 100% in the estimate for

*for ocean entry year 2007, even before considering the sampling variance for*

*R*_{2.2}/_{2}*S*

_{2}(which was ± 50%). We concluded that the 2007 data for age-2 smolts and their returns provided a poor estimate of smolt survival that could not be compared with the unusually low survival of age-1 smolts that also entered the ocean in 2007.

## TABLE 3.

Age-2 Sockeye Salmon smolt composition in returns, 2000–2009. The 95% binomial confidence limits (BCLs) are based on age counts considered separately within each age sample, as opposed to a multinomial confidence limit. *R*_{2}/*R*_{1} is the ratio of returns from age-2 smolts (*R*_{2.1} + *R*_{2.3}) to returns from age-1 smolts (*R*_{1.1} + *R*
_{1.2} + *R*_{1.3}).

### Comparison of *SAR*
_{1} and *SAR*
_{2}

The accepted measure of Chilko Lake Sockeye Salmon smolt survival, *SAR* _{1}, shows substantial year-to-year variability as well as decadal trends (Figure 3A). The mean of age-2 smolts to adult returns (*SAR*
_{2}) was 16% higher than for *SAR* _{1} (Table 1), but this was not statistically different (paired r-test: *P* = 0.20; when log transformed, *P* = 0.57). Neither *SAR*
_{1} nor *SAR*
_{2} was correlated by brood year (n = 45, df = 35, * r^{2}
* = 4%,

*P*= 0.24), but both were weakly correlated by ocean entry year (n = 45, df = 39,

*r*

^{2}= 18%,

*P*= 0.09). Deleting the 2007 outlier in

*SA*had little effect. Both smolt survival time series (Figure 4A) showed similar trends: survival generally rose from 1960 until 1990 and decreased thereafter. Based on the determination that one GAM smoother was as good as two (ANOVA:

*R*_{1}*P*[>F] = 0.29), the low-frequency (decadal) signals in

*SAR*

_{1}and

*SA*were indistinguishable. The temporal trends for

*R*_{2}*SAR*were calculated from log-transformed data (Figure 4A, B) because those GAM models had normally distributed residuals. Although

*SAR*

_{1}and

*SA*were weakly correlated, the first principal component explained 71% of the variance of both despite some of the

*R*_{2}*SA*estimates being poor (but note that when

*R*_{2}*r*

^{2}= 0, PC

_{1}will capture 50% of the variance). These results encouraged us to find an estimate for

*SAR*that combined the estimates of

*SAR*

_{1}and

*SA*, as described in the Appendix.

*R*_{2}###
**Weighting SA***R*_{2} and combining SA*R*_{1} and SA*R*_{2}

*R*

_{2}and combining SA

*R*

_{1}and SA

*R*

_{2}

The age distribution of smolts by ocean entry year (*OEY*) and the age distribution of their survivors were compared as logits, plotting l*og( R_{2.2}/R*

_{1.2}) against log(

*S*

_{2}/

*S*

_{1}) The result (Figure 5 A) confirmed our hypothesis that the age distribution of smolts was consistent with the age distribution of returns for those years when

*S*

_{2}and their survivors,

*had relatively low sampling variance. When*

*R*_{2.2},*S*

_{2}was below 3% of

*S*

_{1}(below —1.5 on the xaxis of Figure 5 A), the scatter began to increase, indicating that the relationship was being lost due to high sampling variance. This result was in agreement with the binomial sampling variance for the cases where the sample sizes were known (Tables 2, 3), e.g., below 20/700 for

*The median weight (Table 4; Appendix) was 0.42, so many cases were strongly downweighted (Figure 5B, D). Low weights for ocean entry years 1965 (related to visual counts in 1967), 1980, 1992, and 2007 essentially removed their influence. The strongest downweighting (six lowest weights) were for*

*R*_{2.2}.*SA*values below 0.05, so one of the effects of weighting was to remove spuriously low values of

*R*_{2}*SA*that may be due to underestimates for

*R*_{2}*R*

_{2.2}. High values of

*SAR*

_{1}in 1969 and 1972 were outliers compared with adjacent years and these were reduced by including

*SA*(Figure 3B). Not all the extreme values of

*R*_{2}*SA*were strongly downweighted; the six highest values of

*R*_{2}*SA*had weights near the median weight. High values of

*R*_{2}*SA*with high weights in 1986 and 1987 were particularly influential and reinforced the pattern of high

*R*_{2}*SAR*

_{1}in the late 1980s. Correlation between

*SA*and

*R*_{2}*SAR*improved after applying these weights (Figure 6) from

_{1}*= 17% to*

*r*^{2}*r*

^{2}= 22% and log transformations increased the weighted correlation to

*r*

^{2}= 33%. Deleting all cases with below-median weight had a similar effect on the correlation (from

*r*

^{2}= 18% to

*r*

^{2}= 22%). Deleting the cases with the lowest 25% of weights increased the proportion of

*SAR*

_{1}and

*SA*variance explained by the first principal component from 71% to 83%.

*R*_{2}We examined the effect on *SAR* of downweighting *SAR* _{1} by 50% for years with low values of *R*_{1.2}as suggested by Bradford et al. (2000). There was almost no effect because *SAR*
_{1} and *SA R_{2}
* estimates were similar in most of those ocean entry years (1963, 1967, 1975, 1979, 1983, 1987; see Table 4) although the value for 1987 increased from 13.0% to 14.6%. We did not weight

*SAR*

_{1}or

*SAR*on this basis, and we had no other indicators for the reliability of

*R*

_{1.2}estimates.

### Density Dependence

The conventional fit of *R*_{1} and *S*_{1} to a Ricker curve (equation 4; Figure 7A) indicated that density dependence was a significant effect in *S*_{1} survival (*r*^{2} = 25%, *P* = 0.0002 when the 2007 outlier is included). The resulting estimates for Ricker parameters a and ß (0.11 and —0.024, respectively) were sensitive to the 2007 outlier, such that maximum returns were 1.7 × 10^{6} at 41 × 10^{6} *S*_{1} before deleting the 2007 case and 2.7 × 10^{6} at 77 × 10^{6} *S*_{1} after. The regression was less convincing after correcting for autocorrelation (*r*^{2} = 7.8%, *P* = 0.070). As expected because of low densities, * R_{2}
* and

*S*

_{2}considered by themselves (Figure 7B) did not show density-dependent survival

*(P*= 0.13). Our novel consideration of

*survival in relation to*

*S*_{2}*S*

_{1}density (Figure 7C) also did not show density-dependent survival

*(P =*0.17). We also determined that the slope from a regression of log

*R*on log

*S*equation (9) was statistically indistinguishable from 1.00, indicating no density dependence (this result is described below).

### Regime Effects

Survivals based on the *SAR* estimates (*SAR* _{1} and *SA R_{2}
* combined) differed among four ecosystem regimes (ANOVA based on equation 7:

*P*< 0.001). Mean

*SAR*(as percent) increased from 8.6% during 1960–1976 to 10.7% during 1977–1988 and 10.2% during 1989–1998, and then decreased to 4.1% during 1999–2008 (Table 5). After correction for multiple comparisons, the only significant pairwise differences indicated that 1999–2011 survival was lower than for 1977–1988

*(P =*0.008) and 1989–1998 (

*P*= 0.028). Adjusting the regimes to change at 1991 instead of 1989 showed that

*SAR*(in percent) dropped by 4.3% (2.2) between 1977–1990 and 1992–1998, but this was not a significant difference

*(P*= 0.17).

Fitting a weighted smoother to the smolt and adult returns time series (equation 8; Figure 4C) indicated a trend of increasing survivals from 1960 to 1990 without a clear change at 1977, then a trend of lower and decreasing survivals from 1992 to 2008 (there is no estimate for 1991). Survival effects from ecosystem regime changes in 1977 and 1999 appeared, in this case, to be artifacts from describing continuous trends as discrete means.

## TABLE 4.

Percent survival from smolt to adult returns by ocean entry year *(OEY)* for age-1 and age-2 Sockeye Salmon smolts (*SAR* _{1} and *SA R_{2}
*), the estimated reliability of

*SA*expressed as a weight ranging from zero to one, and

*R*_{2}*SAR*calculated as the weighted geometric mean of

*SAR*

_{1}and

*SAR*

_{2}.

### Fisheries Management

As well as the change in *SAR* after *OEY* 1991, the effects of fisheries management created large changes to postlake survival of Sockeye Salmon beginning in the 1990s (Table 6). Until 1990 there were consistently high losses between recruits and spawners (Ψ) and a consistent 4-year cycle of low returns and low escapements (1953, 1957,… 1985, 1989) (Figure 2E, F). Starting in 1990, escapement suddenly increased and subsequently remained high for every year with the exception of 2004. The mean escapement from 1990 to 2008 (0.33 × 10^{6}) was 250% of that from 1960 to 1989 (0.13 × 10^{6}). There were no low escapements comparable with the 1953 cycle line for 14 years after 1989. Starting in 1995, Ψ decreased from about 80% to 50% and remained low until 2009, our last estimate. In 2009, approximately 94% of adult returns from the ocean reached Chilko Lake because the fishery was immediately closed when the record-low smolt survival from ocean entry year 2007 was recognized (Figure 2F).

## TABLE 5.

Mean and SE for *SAR* (%) of Sockeye Salmon for ocean climate regimes beginning in 1977, 1989, and 1999. The last two rows are similar to rows 2 and 3, but assume that an ocean survival regime started in 1991 instead of 1989, as indicated by the trends in *SAR.*

### Smolt Length, Condition, and Age

Mean lengths and weights of age-1 smolts (*S*_{1}) were not significantly different among time periods (Table 7). On average, *S*_{2} were longer and correspondingly heavier after 1991, largely due to *S*_{2} in 2007, when the mean length of *S*_{2} was 160 mm. During the 1977–1991 period of high survivals, *S*_{1} were in significantly better condition than during the 1992–2009 period of low survivals *(t-*test: *P* = 0.003). Condition of *S*
_{2} was also better before 1991 than after *(P* = 0.06).

The "base case" for analyzing the effect on smolt survival of smolt quality (age, length, and condition) and time trends was the plot of log_{e}(R) versus log_{e}(S) corresponding to equation (10), Figure 8, and model A in Table 8. The 2007 outlier was deleted for this analysis. Plots of partial autocorrelation for the residuals of the full model (row D in Table 8) showed that autocorrelation was not significant, and fitting a first-order autocorrelation model resulted in a small value (Φ = 0.024) for off-diagonal elements of the residuals covariance matrix. We proceeded without including a correction for autocorrelation in these regressions. Heterogeneous sampling variance in * R_{2}
* and

*S*

_{2}was explored further by using our weights (Figure 5D) as a variance covariate (a predictor of the variance of residuals but not a predictor of

*Zuur et al. 2009) to fit new estimates of sampling variance in a generalized least-squares regression. The result was similar to our weights, but suggested a smaller variance (higher weights) for the most reliable cases of*

*R*_{2}\*and*

*R*_{2}*S*

_{2}. This supported our parameterization of sampling variance, and we proceeded with our weights as the more conservative approach. Fitting autocorrelation and heterogeneous variance models simultaneously showed little interaction. Plots of residuals (adjusted for weights) against fitted values for the full model showed no trends, and the residuals had a near-normal distribution.

## TABLE 6.

Prespawning loss (primarily exploitation), Ψ, and female escapement statistics for Sockeye Salmon before and after 1994.

## TABLE 7.

Length (mm), weight (g), and condition (%) of Sockeye Salmon smolts by age and by selected ranges of ocean entry years. Weights are from means of preserved samples. Condition is deviation (%) of observed from predicted mean weight based on means of preserved lengths.

Covariates such as age and length were interpreted as a multiplicative effect on survival rates in these models. We expected that accounting for the effects of climate trends and smolt abundance would reveal patterns in the residual variance that were due to smolt quality; however, we could not discern a significant effect from smolt age, length, or condition in the regressions with various combinations of covariates (Table 8). Linear temporal trends broken at 1991, without effects from smolt quality (trial B in Table 8), provided the best model as identified by the AIC and explained 91% the returns across both age-classes of smolts (Figure 4D). For this model, the fitted parameter for log_{e}(S) was 0.955 (0.033), within two SEs of 1.00, implying that there was no curvature in the smolts-to-adult returns relationship. When the parameter for log_{e}(S) was forced to 1.00 (trial G), the change in *r*^{2} and AIC was small and the difference between 0.95 and 1.00 for this parameter could have arisen due to sampling variance (ANOVA between trials B and G: *P =* 0.17).

## TABLE 8.

Regression results from examining smolt age, smolt length, and time trends as predictors of Sockeye Salmon smolt survival. In the model notation, *S* is log_{e} (smolts) of both ages, *R* is Iog_{e} (returns) from those smolts, *L* is the mean length of smolts (by age), A is a factor for smolt age (1 or 2), *T* is a factor indicating whether ocean entry year *(OEY)* is before or after 1991 for fitting two trends with time, and exp(S) is the Ricker term for density dependence. The number of cases is 93 for all models, weights for * R_{2}
* and

*were applied, and the 2007 outlier for*

*S*_{2}*R*

_{1}was excluded. Trial D is the full model that was used to examine autocorrelation and other patterns in residuals. Trials with smolt condition (K) were deleted from this table because it was never a significant effect when added to these models.

## TABLE 9.

Abundance of Sockeye Salmon returns by age-groups (mean, SD, CV [SD/mean]) and correlations between returns within a brood but after differing lengths of time at sea. Degrees of freedom (df) for correlations reflect correction of probabilities for autocorrelation.

The model in trial G predicted log_{e}(SAR) without density dependence and the parameter estimates (with SEs) were log_{e}(R/S) = -58.81(14.69) + 0.02857 (0.00744) × *OEY* before 1991, and l*og _{e}(R/S)* = 92.59 (43.21) – 0.04777 (0.0217) ×

*OEY*after 1991. When we added — ßS, the Ricker term for density dependence (trial H), the change in

*and AIC was small and the difference was not significant (ANOVA between trials G and H:*

*r*^{2}*P*= 0.21). After removing temporal trends, intrastock density-dependent survival of Chilko Lake Sockeye Salmon smolts was undetectable.

### Marine Survival after the Ocean Entry Year

Precocious males (i.e., jacks, *R*_{1.1} or *R*_{2.1}) were poor predictors of the subsequent abundance of their siblings, * R_{1.2}
* and

*or*

*R*_{1.3}*and*

*R*_{2.2}*(Table 9). There was a significant positive correlation between*

*R*_{2.3}*R*

_{1.2}and subsequent

*R*

_{1.3}, both of which were abundant and estimated relatively precisely. Similarly,

*and*

*R*_{2.2}*although less abundant, were also correlated. We pursued the correlation between*

*R*_{2.3},*and*

*R*_{1.2}*, noting that it was affected by three factors: (1) shared dependence on*

*R*_{1.3}*S*

_{1}, (2) year-to-year variation in the proportions that return after 1, 2, and 3 years at sea (a separate effect from survival), and (3) shared mortality in the first 2 years after leaving Chilko Lake. After removing the contribution of

*S*

_{1}to both

_{1.2}and

*R*

_{1.3}, both time series shifted to lower survivals after 1991, but numerous independent outliers in both time series resulted in a nonsignificant correlation (P = 0.29). There was no evidence of an effect from a common ocean entry year between successive returns within broods.

### Changes in Sea Residence

We examined how the age distributions of returns varied with time (Figure 9) because year-to-year variability in duration of sea life contributed to the residual variance of the regressions described above. Duration of marine life for age-1 smolts increased from 2.015 (0.005) years during 1948–1967 to 2.092 (0.007) years during 1989–2008, a 3.8% (0.88) change in 40.5 years. This effect was almost entirely because the proportion of *R _{1.3}
* increased 350% from 2.7% (0.84) in 1948–1967 to 9.5% (1.16) in 1989–2008. There was a linear trend of 30.5% (10.0) per century in the proportion of

*R*

_{1.3}from 1974 to 2009 (Figure 9B).

In the same period, the proportion *of R*
_{1.1} (jacks) decreased by 75% from 1.2% (0.21) to 0.3% (0.05) (Figure 9A). The equivalent pattern of older returns from *S*_{2} is less clear (Figure 9D); however, the proportion of returns from *S*_{2} that were *R*_{2.1} (jacks) dropped to nearly zero after 1994, whereas proportions greater than 10% had been frequent previously (Figure 9C).

### The Marine Environmental Signal

Corresponding to the temporal pattern of *SAR* for Chilko Lake Sockeye Salmon (Figure 4D), there was a break in trends near 1991 for many indices of the density of salmon in the North Pacific Ocean that might have been competing with Chilko Sockeye Salmon (columns 1 [North America] and 3 [Asia] of Figures 10, 11). For North America, this included catches of Sockeye and Chum Salmon, and hatchery releases of Sockeye, Pink, and Chum Salmon, but Pink Salmon catches continued to increase after 1991. For Asia, this included catches of all three species, and the trends broke near 1987 for Sockeye Salmon and 1994 for Chum Salmon. Asian hatchery releases of Sockeye Salmon and particularly the massive releases of Chum Salmon changed near 1991, but Pink Salmon releases were variable and did not show a clear break. Reliable estimates of hatchery releases of Sockeye Salmon by Russia (and therefore Asia) before 1991 were not available, but the numbers released were probably low. Climate indices (columns 1 and 3 in Figure 12), in contrast to salmon density indices, did not show a pattern break near 1991 with the exception of PDO (winter and summer) and ENSO. Coastal indices of upwelling, SST, and Strait of Georgia Pacific Herring did not have trends with a break near 1991.

To compare trends in *SAR*, with its substantial year-to-year variance, with trends in climate and salmon density, we examined the correlation of *SAR* with each of the 24 indices, before and after 1991 (columns 2 and 4 in Figures 10–12). For North America, there were weak but significant correlations (at α = 0.05) with catches (column 2 in Figure 10) and hatchery releases (column 2 in Figure 11) of competing salmon species before 1991 (except for catch of Chum Salmon and hatchery release of Sockeye Salmon) but only Chum Salmon hatchery releases after 1991 were close to significance (P = 0.07). For Asia there was a significant correlation of l*og _{e}(SAR)* with catches of Sockeye Salmon before 1991 (column 4 in Figure 10) that suggested interstock density dependence, but evidence of this after 1991 was weak. Asian Chum Salmon catches were also correlated with lo

*g*before 1991 but the slope was positive, indicating an effect on both from increasing ocean productivity, but that relationship disappeared after 1991. Sockeye Salmon catches in Asia before the early 1980s were primarily of Russian and North American fish by Japanese fishers, but afterwards, Russia caught increasing proportions of the catch, and now virtually all the Sockeye Salmon caught in Asia are by Russia (Irvine et al. 2012). The only evidence of

_{e}(SAR)*SAR*being depressed after 1991 by hatchery releases in Asia was a marginal correlation with Pink Salmon (column 4, row 2 in Figure 11; P = 0.053). The only significant correlation of

*SAR*with climate and coastal indices that we examined was from coastal upwelling during summer at 48°N after 1991, despite there being no trend in that index and no correlation before 1991 (column 2, row 4 in Figure 12).

## DISCUSSION

Results from our examination of 50 years of information on two age-classes of Sockeye Salmon smolts emigrating from and returning to Chilko Lake helped us understand reasons underlying variability in abundance and survival for Sockeye Salmon in general. When arranged by ocean entry year, survivals of age-1 and age-2 smolts were weakly positively correlated, suggesting that each was a noisy estimate of survival. To understand the error associated with these survival estimates, we used a novel principal components approach to estimate the relative sampling variance for uncommon age-2 smolts and their survivors, allowing a weighted combination of survival estimates from the two smolt ages. The combined time series reduced the sampling variance in the original age-1 smolt survivals sufficiently to indicate a linear trend of increasing smolt survival during 1960–1990 that suddenly changed at or near 1991 to a lower and declining trend in survival during 1992–2008.

The judicial inquiry on Fraser River Sockeye Salmon was precipitated in part by exceptionally low Sockeye Salmon returns in 2009 (Cohen 2012a, 2012b, 2012c). The survival of Chilko Lake Sockeye Salmon that went to sea as a record high abundance of age-1 smolts in 2007 and returned as *R*_{1.2} in 2009 was anomalously low in our analyses (i.e., an outlier), consistent with the extremely low returns that year for Fraser River Sockeye Salmon in general. We were unable to compare this low survival of age-1 smolts with the survival of age-2 Chilko Lake smolts that also entered the ocean in 2007 and returned in 2009 because of the sampling variance associated with low counts for age-2 smolts within the age samples. In contrast to Fraser River Sockeye Salmon, survivals for many non-Fraser River Sockeye Salmon (e.g., Columbia River, west coast of Vancouver Island, central and northern British Columbia) that returned in 2009 were not abnormally low (K. D. Hyatt, DFO, Pacific Biological Station, personal communication). It appears that the low survivals of Fraser River Sockeye Salmon returning in 2009 were primarily a consequence of 2007 marine conditions in the Strait of Georgia that also affected other species including Pacific Herring (Beamish et al. 2012), as well as conditions to the north of the Strait of Georgia, especially within the Queen Charlotte Sound-Hecate Strait region (Thomson et al. 2012). The high returns to Chilko Lake in 2010, coincident with the highest returns to the Fraser River watershed in a century, were the consequence of high egg-to-smolt survival in freshwater yielding large numbers of smolts going to sea in 2008 followed by average (compared with long term) smolt-to-adult survivals.

Insignificant differences between survivals of age-1 smolts and the larger age-2 smolts were unexpected. Koenings et al. (1993) found that 30% of the variation in Sockeye Salmon smolt survival was explained by smolt size, much more than smolt age, as well as a south-to-north cline of increasing survivals. When Henderson and Cass (1991) compared scale-based growth rates for 3 years of age-1 Chilko Lake Sockeye Salmon smolts to those of returning adults, there was evidence that larger age-1 smolts survived better than smaller smolts in the same cohort. However, as in our study, Henderson and Cass (1991) found no significant effect of mean size on the survival of age-1 smolts, based on 34 years of data. In a more recent examination of growth patterns for Fraser River Sockeye Salmon, McKinnell and Reichardt (2012) found that early marine growth, measured as size at the end of the first marine growing season of returning adults, had a minor effect on interannual survival differences. It appears the complex relationship between smolt size and survival is affected by many factors including smolt age and timing and stream latitude.

Our finding of increased sea residence for Chilko Lake Sockeye Salmon agreed with the findings of others who have documented increasing numbers of older Sockeye Salmon recruiting in recent years (e.g., Holt and Peterman 2004). Increased sea residence was consistent with a reduction in first-year marine growth for Fraser River Sockeye Salmon commencing in 1977, which occurred in spite of concurrent increases in biological productivity in the coastal marine zone (McKinnell and Reichardt 2012).

A long-standing puzzle in Fraser River Sockeye Salmon biology is the explanation for cyclic dominance, where some Fraser River Sockeye Salmon populations exhibit a pattern of strong returns every 4 years, often with a subdominant year and two weak years in between, while other populations do not (e.g., Ricker 1950; Welch and Noakes 1990). In the case of Bowron Lake (mid-Fraser River watershed), Sockeye Salmon demonstrated a pattern of cyclic dominance from 1959 to 1982, but not before or after that time; reasons for this shift are not clear (Walters et al. 2004; Grant et al. 2011). We found that before 1990, Chilko Lake Sockeye Salmon exhibited clear cycles marked by a weak return every 4 years (i.e., 1961, 1965, ... 1989). This pattern changed suddenly after 1989 when spawner escapement doubled as a result of unusually high recruitment during 1990– 1994 (predominately *R* _{1}, Figure 2C) and then reduced fishing during 1995–2009 (Ψ, Figure 2E). The pattern of cyclic dominance before 1990 influenced our weights for survival estimates for age-2 smolts. During the period of cyclic dominance, large escapements one year followed by small escapements the next year resulted in a large proportion of age-2 smolts comigrating with relatively few numbers of age-1 smolts. This age distribution was repeated in the adult returns 2 years later. We used that consistency in age distributions to parameterize sampling variance for *SA R_{2}.* In ocean entry years with a high ratio of age-2 to age-1 smolts, survival estimates for age-2 smolts were the most reliable and had a strong influence (heaviest weights) on the combined survival estimates

*(SAR).*These weights may be revised if or when additional sample size information becomes available.

Declining survivals for Chilko Lake Sockeye Salmon smolts after 1991 corresponded to decreasing productivity (returns per spawner) noted about that time for many North American Sockeye Salmon stocks outside of central and western Alaska (Peterman and Dorner 2012). We have no real explanation for differences among populations in the timing of declines although stock-specific migration patterns (Tucker et al. 2009) or differences in relative strengths of particular brood years may be important. Our results suggest that this general pattern of decline for Sockeye Salmon resulted from factors downstream from rearing lakes, which is in agreement with findings from a 2010 workshop on Fraser River Sockeye Salmon (Peterman et al. 2010). Although the influence of factors between the lake and the ocean on survival estimates cannot be excluded, the wide geographic distribution of Sockeye Salmon populations with similar patterns in productivity (Peterman and Dorner 2012) suggests that freshwater factors are unlikely to have played a major role in the widespread decline.

The onset of declining survival about 1991 also corresponded approximately to a shift in climate regimes from a period of above-average productivity for many populations of Pacific salmon to one that was less productive (e.g., Hare and Mantua 2000; Irvine and Fukuwaka 2011). Our analysis did not suggest breaks in the survival time series corresponding to ecosystem regime shifts in 1977 and 1999.

After removing temporal trends, we found no evidence of an intrastock density-dependent effect on smolt-to-adult survival. Similarly, McKinnell and Reichardt (2012) found no relationship between the numbers of smolts leaving Chilko Lake and their subsequent first year of marine growth. In the ocean, individuals from any one salmon population will potentially compete with fish from many other populations, so one might expect density-dependent effects in the ocean to be due to interstock competition. A logical extension of this concept would be to replace estimates of Chilko Sockeye Salmon smolt numbers with estimates of Sockeye Salmon smolts from the Fraser River and perhaps other rivers draining into the northern California Current. In this paper we extended the concept even further to include interspecies competition and investigated correlations between Sockeye Salmon survival and abundance estimates of Sockeye, Pink, and Chum Salmon in the eastern and western North Pacific Ocean.

It is interesting that trends in Asian Sockeye Salmon catches reversed in the 1980s (row 1, column 3 in Figure 10) and that these trends were correlated with smolt survivals in Chilko Sockeye Salmon (row 1, column 4 in Figure 10). Declining catches before 1990 were chiefly the result of high catches of Russian and North American Sockeye Salmon in the high seas by Japan early in the time series. In contrast, increasing catches after 1990 were primarily the result of increasing catches of Russian Sockeye Salmon by Russia. Are Russian fisheries the real reason for the apparent declines in survivals of North American Sockeye Salmon after 1990? Using data from Irvine et al. (2012, their Tables 1 and 7), we determined that since 1990, Sockeye Salmon catches by Russia have been increasing and in recent years have constituted 25–30% of the total North Pacific Sockeye Salmon catch. These estimates may be low; Dronova and Spiridonov (2008) reported that Sockeye Salmon catches in Kamchatka were 1.5–3 times higher than reported. However, Russian salmon fisheries primarily occur close to Russia, and although some Canadian Sockeye Salmon are found in the Bering Sea (Habicht et al. 2010; Beacham et al. 2011), their numbers in the western North Pacific Ocean are low. We therefore think it is unlikely that declining North American Sockeye Salmon survivals after 1991 are a consequence of these salmon being caught in Russian fisheries and not included in marine harvest estimates. As an alternate explanation, the correspondence of Chilko Lake Sockeye Salmon survival to the patterns of North American and Asian catches and hatchery releases suggests that density-dependent processes in the North Pacific may be affecting Sockeye Salmon survival.

Peterman and Dorner (2012) suggested that mechanisms controlling Sockeye Salmon marine survival patterns would likely (1) operate at large spatial scales or where many Sockeye Salmon populations overlap; (2) affect populations from Puget Sound to southeastern Alaska similarly, while having an inverse effect on stocks from central and western Alaskan; or (3) have been present historically, but intensified recently. With this in mind, and using results from our examination of correlations between various data sets and our new time series supplemented by results from the literature, we developed a conceptual model to explain how changes in ocean climate and salmon species abundance might control the survival and growth of Sockeye Salmon (Table 10).

In our model, we suggest that increasing biological productivity before 1991 resulted in greater carrying capacity for salmon. Also during this period, expanding hatchery programs yielded greater Pink, Chum, and Sockeye Salmon biomass. Responses to more salmon competing for more resources included generally improved survivals for Sockeye Salmon but negative responses in terms of growth (Table 10). Survival was essentially independent of growth during this period. After 1991, reduced ocean productivity resulted in declining carrying capacity. Continuing increases in Pacific salmon numbers and biomass, largely from hatcheries, further intensified inter- and intraspecific competition for food. In response, Sockeye Salmon growth rates continued to decline, and as a result of many fish not achieving some critical size, survivals began to decline.

Evidence to support this model comes from various sources. Kaeriyama et al. (2009) concluded that Pacific salmon carrying capacities in the North Pacific Ocean peaked between 1985 and 1994 for Sockeye, Pink, and Chum Salmon, presumably alleviating competition effects until, in the case of Chilko Lake Sockeye Salmon, approximately 1991. We showed earlier that release numbers from salmon hatcheries generally increased until about 1991 (columns 1 and 3 in Figure 11), and the literature provides widespread evidence of increasing Pacific salmon abundance and survivals between 1977 and 1991 (e.g., Beamish and Bouillon 1993; Mantua et al. 1997). With respect to the three primary salmon species, Chum Salmon consistently constitute the largest biomass in the North Pacific (smaller Pink Salmon are the most numerous), followed by Pink Salmon and then Sockeye Salmon (Eggers 2009). Japanese Chum Salmon migrate seasonally between the Bering Sea (summer-fall) and the Gulf of Alaska (winter-spring) (Urawa et al. 2009). Sockeye Salmon from Puget Sound to southeastern Alaska spend most of their time within the Gulf of Alaska where they overlap with, and presumably compete with many species, including Asian Chum Salmon during the winter. Recent genetic evidence also demonstrates North American Sockeye Salmon extend farther northward and westward into the Bering Sea than thought previously where they overlap with Asian Sockeye and Pink Salmon (Habicht et al. 2010). Body size for many Sockeye Salmon populations increased from the 1960s until the late 1970s, followed by decreases (Eggers and Irvine 2007). Similarly, McKinnell and Reichardt (2012) reported a decline in early marine growth from 1970 to 1990 of Chilko Lake and Birkenhead River Sockeye Salmon, based on scale growth of survivors. We suggest that before about 1991 expanding numbers of Pacific salmon were increasingly competing for a healthy but still finite food resource, and this view is consistent with the overall finding of reduced growth rates but not reduced survival for Sockeye Salmon (Table 10).

## TABLE 10.

Conceptual model illustrating possible impacts of changing marine productivity, hatchery releases, and hatchery effectiveness on survival and growth responses of Sockeye Salmon from British Columbia and southern Alaska. Up arrows “↑” and down arrows “↓”indicate increases and decreases while diagonal arrows “↗” and right arrows“→” indicate minor increases and little change, respectively.

Kaeriyama et al. (2009) concluded that carrying capacities for Pacific salmon declined after the mid-1990s. Yet continued increases in hatchery salmon releases in Asia after 1991, as well as improved survival for many hatchery fish resulting from the release of larger or better-adapted young salmon (Dushkina 1994; Morita et al. 2006; Kaev and Ignatiev 2007), resulted in increasing salmon biomass. Proportions of hatchery-origin fish in the North Pacific have been increasing since 1990 and constitute 50–62% of the total Chum Salmon, 10–13% of the Pink Salmon, and 4–10% of the Sockeye Salmon in the North Pacific Ocean (Ruggerone et al. 2010; Kaeriyama et al. 2012). Eggers (2009) estimated that at least 39% of the salmon biomass in a recent 10-year period was made up of hatchery-origin Pink and Chum Salmon. We suggest that after 1991 more salmon were competing for a reduced food supply, resulting in greater competition than before 1991 (Table 10).

Chum and Pink Salmon, which constitute the bulk of the salmon in the North Pacific Ocean, appear to have ecological advantages over Sockeye Salmon during periods of increased competition. For instance, Kaeriyama et al. (2012) found Chum Salmon were better able to alter their diet during periods of changing prey communities than were other salmon species. Based on a detailed review of the literature, Ruggerone and Nielsen (2004) suggested that Pink Salmon are the dominant competitor among the salmon species, affecting other species by reducing the availability of prey. Ruggerone et al. (2003) found reduced Sockeye Salmon growth during their second and third year at sea in years when Asian Pink Salmon were most abundant.

After 1991 there was no evidence of further declines in Fraser River Sockeye Salmon growth rates, although the early marine growth of fish surviving to return was less than average due to earlier declines (McKinnell and Reichardt 2012). Perhaps growth rates continued to decline after 1991, but slow-growing Sockeye Salmon died at higher rates than did faster-growing Sockeye Salmon (Ruggerone et al. 2005; Farley et al. 2011 ) such that measured growth rates in the survivors did not decline. We suggest that Sockeye Salmon survivals began to decline after 1991 for many stocks because an increasing portion of the juveniles were unable to grow to the size (or possess sufficient energy reserves) required to survive some life history stanza, similar to the critical size hypothesis (e.g., Beamish et al. 2004a). Further evidence of increasing competition is the increased proportion of Chilko Lake Sockeye Salmon returning late after three winters at sea instead of two. This linear trend started in the 1970s when 2% returned late, and by 2009, 12.6% returned late (Figure 9B), presumably as a result of slower growth.

Our results do little to clarify which life history stage is being affected by density-dependent effects in the ocean. The first year of marine growth for Chilko Sockeye Salmon is greater during even years, when young Pink Salmon are abundant, than in odd years when they are not (McKinnell and Reichardt 2012), but it is not known whether this pattern continues during the Sockeye Salmon's life. Correlations between survival and abundance (catch) of Pink Salmon, lagged to align with ocean entry year, might be interpreted to indicate that density-dependent effects occurred in the first 18 months at sea for Chilko Sockeye Salmon, the marine residence period for Pink Salmon. However, older Sockeye Salmon could also compete with Pink Salmon. There were also correlations between Sockeye Salmon survival and abundances of Sockeye and Chum Salmon as well as hatchery releases, and these fish would have shared the ecosystem with Chilko Sockeye Salmon throughout their lives. Examining correlations at lags of 1, 2, and 3 years may help understand this issue. Ruggerone et al. (2005) hypothesized that large numbers of Pink Salmon depressed zooplankton and other prey during odd-numbered years, thereby reducing the growth rates of Bristol Bay Sockeye Salmon, primarily in their second ocean year.

Our conceptual model (Table 10) satisfies the first (large spatial scales) and third (present historically but intensified recently) requirements of Peterman and Dorner (2012). Although the linkages proposed here among ocean productivity, hatchery releases, and Sockeye Salmon survival and growth are simplistic, this model offers promise for a more complete explanation of how salmon population dynamics are linked to North Pacific ecosystem dynamics. Further research, including the spatial and dietary overlap of marine salmon and consideration of nonsalmon competitors, is required to test the hypotheses of this model.

In conclusion, while both local- or regional-scale and oceanbasin-scale processes are important in controlling survival patterns for Pacific salmon, our detailed examination of data from one population suggests that density-dependent processes operating far from coastal North America may be important in determining long-term survival patterns for Sockeye Salmon. The importance of regional effects was evidenced by the anomalously low survival of smolts from ocean entry year 2007. We support the approach advocated by Peterman and Dorner (2012) of simultaneously evaluating multiple salmon populations, but also stress the importance of detailed examinations of populations such as Chilko Sockeye Salmon where mortality can be partitioned between events occurring primarily in the ocean versus events that occur in freshwater. We encourage other researchers to closely examine their data sets for additional time series of information that may be available, as exemplified by the relatively uncommon older smolts that we were able to include in this analysis.

## ACKNOWLEDGMENTS

Information on Chilko Lake Sockeye Salmon result from decades of collaborative work by staffkk of Fisheries and Oceans Canada (DFO) and the Pacific Salmon Commission (PSC). We thank all those who have been involved in the collection of data on Chilko Lake Sockeye Salmon, and especially Sue Grant (DFO) and Mike Lapointe (PSC) who provided data and comments on an early manuscript. Discussions with David Welch of Kintama Research led to improvements in text and analysis. We appreciate data provided to us on Strait of Georgia Pacific Herring by Jake Schweigert (DFO) and on Fraser River Sockeye Salmon by Keri Benner, Tracy Cone, Jeremy Hume, and Sue Grant (DFO) and Mike Lapointe and Steve Latham (PSC); comments on the manuscript by Brendan Connors; and the editorial assistance of Lana Fitzpatrick (DFO). The paper benefited from thoughtful and constructive comments by two anonymous reviewers.

## REFERENCES

*in*C. C. Krueger and C. E. Zimmerman , editors. Pacific salmon: ecology and management of western Alaska's populations. American Fisheries Society, Symposium 70, Bethesda, Maryland. Google Scholar

*Oncorhynchus nerka*) wild salmon policy status using abundance and trends in abundance metrics. Canadian Science Advisory Secretariat Research Document 2011/087. Google Scholar

*Oncorhynchus nerka*). Canadian Journal of Fisheries and Aquatic Sciences 48:988–994. Google Scholar

*Oncorhynchus nerka*) in a changing environment. Canadian Journal of Fisheries and Aquatic Sciences 61:2455– 2470. Google Scholar

*Oncorhynchus nerka*)

*:*effects of smolt length and geographic latitude when entering the sea. Canadian Journal of Fisheries and Aquatic Sciences 50:600–611. Google Scholar

*Oncorhynchus nerka*) in relation to juvenile Pink (

*Oncorhynchus gorbuscha*) and Sockeye salmon abundance. Canadian Journal of Fisheries and Aquatic Sciences 69:1499–1512. Google Scholar

*Oncorhynchus nerka*). Canadian Journal of Fisheries and Aquatic Sciences 41:1825–1829. Google Scholar

*Oncorhynchus nerka*) populations in western North America. Canadian Journal of Fisheries and Aquatic Sciences 69:1255–1260. Google Scholar

*Oncorhynchus nerka*) stocks. Canadian Journal of Fisheries and Aquatic Sciences 55:2503–2517. Google Scholar

*Oncorhynchus nerka.*PLoS (Public Library of Science) ONE [online serial] 7(4):e34065. Google Scholar

*Oncorhynchus nerka*) in relation to competition with Asian Pink Salmon (

*O. gorbuscha*) and the 1977 ocean regime shift. U.S. National Marine Fisheries Service Fishery Bulletin 103:355–370. Google Scholar

*Oncorhynchus gorbuscha*) over other salmonids in the North Pacific Ocean. Reviews in Fish Biology and Fisheries 14:371–390. Google Scholar

*Oncorhynchus*gorbuscha) and Alaskan Sockeye Salmon (

*O. nerka*) in the North Pacific Ocean. Fisheries Oceanography 12:209–219. Google Scholar

*Oncorhynchus nerka*) escapement. Canadian Manuscript Report of Fisheries and Aquatic Sciences 2428. Google Scholar

*Oncorhynchus nerka*). Canadian Journal of Fisheries and Aquatic Sciences 47:838–849. Google Scholar

## Appendices

### Appendix: A Proxy Model for Variance of Estimates of *SA**R*_{2}

*R*

_{2}

If two time series have some desired signal in common, but a low correlation because sampling error or other noise has been added to both of them, then they can be usefully combined because the noise will cancel out, but only if the noise is independent. Both *SAR* _{1} and *SA R_{2}
* are unlikely to have completely independent noise because they are from the same field sampling operations; however, they are from different brood years (different parents), have a different freshwater history, enter the ocean at different sizes, and are tracked separately in fisheries statistics. Both time series showed substantial year-to-year variance and little covariance, so an attempt at noise reduction seemed worth pursuing.

We wished to weight estimates of *SA R_{2}
* in inverse proportion to their variance for regressions and statistical tests and wanted to identify cases where either or both of

*S*

_{2}and

*R*

_{2.2}had high sampling variance such that the ratio

*was an unreliable estimate of*

*R*_{2.2}/*S*_{2}*SA*The age distribution of smolts in an ocean entry year, t, should be consistent with the age distribution of returns from those smolts 2 years later, i.e.,

*R*_{2}.from which

and *SAR*_{2} = *SAR*_{1}.

Consistency of age distributions is therefore a prerequisite for survival of age-2 smolts to be comparable with survival of age-1 smolts. When the binomial sampling variance of *S*_{2,t} is large because *S*_{2,t} are a small fraction of the smolts sampled for aging, the associated * R_{2.2,t+ 2}* will also be rare and observed with large sampling variance. Consequently the observed age distribution will be highly variable and the ratio

*will be unreliable.*

*R*_{2.2,t+ 2}/*S*_{2,t}= SA*R*_{2}Our procedure was to plot log_{10} (*R*
_{2.2,t+ 2}/*R*_{1.2,t+2}) versus log _{10} (*S*_{2,t}/*S*_{1,t}) to see whether there was consistency in the age distributions (Figure 5A). This was not a regression situation where we were attempting to predict log _{10}(*R*
_{2}/*R*_{1}) or log_{10}(*S*_{2}/*S*_{1}); rather we were trying to observe the sampling variance in both. We used principal components to describe this relationship. The first principal component, PC_{1}, will be a line describing the consistent relationship in age composition between smolts and returns, if it exists, and will have slope = 1 on this plot. This corresponds to equation (A.2) being true. The second component, PC_{2}, will then capture the scatter of points orthogonal to this line, identifying cases that differ from the hypothesized consistent relationship. Inconsistent cases will have high loadings on PC_{2} and low loadings on PC_{1}. Note that PC_{1} is not simply log_{10}(*R*_{2}/R _{1} ) or *log*
_{10}(*S*
_{2}/S _{1} )• Also, PC_{1} can be used directly as regression weights, *w,* which was how we used it (Figure 5B), and can also be used as a “variance covariate” as a predictor for heterogeneity in the variance of residuals in generalized linear models (e.g., the function gls in the R package nlme as described by Zuur et al. 2009:chapter 4; see also Gelman et al. 2003:chapter 14.6). We explored the possibility of improving our weights by using PC_{1} loadings as a variance covariate for the varConstPower function in GLS models and determined that only a small adjustment resulted.

That this consistency exists is a hypothesis to be tested. Subsequently, examining that consistency on a case-by-case basis will indicate the relative reliability (weight) for each estimate of *SA R_{2}.* Changes in the duration of sea life will interfere in this scheme. We examined the variation in years at sea to see whether this is a problem (Figure 9), but note that sampling variance for age also applies to duration of sea life.

We also examined weights for *SA R_{2}
* as a theoretical variance for the ratio of binomial variables. Assuming that binomial sampling error for the age distribution is the predominant error for

*S*

_{2}and

*, we applied an approximation for the variance of a ratio, derived from the second-order Taylor series approximation (Stuart and Ord 1987). The resulting general formula is*

*R*_{2.2}where *E(S)* is the expectation of the random variable S. In the case of binomial distributions for estimates of age frequency each year, with parameters *p*, *n*, and *q = 1 — p,* the variance of the observed number of *S*_{2} in the age sample of n smolts is

in a specific year (not the variance of the proportion, *pq/n)* and *E(S) =* n*p* is the number of 52 observed in the age sample of *n* smolts. Applying these binomial parameters in equation (A.3), then

Equation (A.3) can be treated more simply than this because var(x) = *E(x)* for the binomial distribution and *pqn* approaches *pn* as *p* approaches 0. So when 52 are a small fraction in the age samples,

The observed counts of 52 appear as squared and cubed terms in denominators in equations (A.5) and (A.6), which is a warning that estimates of *SA R_{2}
* are sensitive to the sampling variance of

*S*

_{2}. We dropped the covariance term in equation (A.3) because it is not measured and thought to be small, despite the above argument that the proportions of

*and*

*R*_{2.2}*S*

_{2}in age samples do covary, so equations (A.5) and (A.6) will overestimate Var

*(SA*As well, we are applying an incomplete model when we estimate

*R*_{2}).*Var*(

*SA*) based on simple binomial sampling; other sources of sampling variance include allocation of catches from mixed stocks of Fraser River Sockeye Salmon to the Chilko Lake stock and application of length-at-age matrices to smolt length distributions. Nevertheless, equation (A.5) was useful as a check for the weights derived from comparing the age distributions of smolts and returns (Figure 5C, D). To apply equation (A.5) when we did not know the size of the age samples, we used the proportions of

*R*_{2}*S*

_{2}and

*R*

_{ 2.2 }in the total smolts and total returns (not the age samples) for

*p,*and we used the average of recent sample sizes for n.

### APPENDIX REFERENCES

Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2003. Bayesian data analysis, 2nd edition. CRC Press, Boca Raton, Florida.

Stuart, A., and J. K. Ord. 1987. Kendall's advanced theory of statistics, volume 1: distribution theory. Charles Griffin and Company, London.

Zuur, A. F., E. N. Ieno, N. J. Walker, A. A. Saveliev, and G. M. Smith. 2009. Mixed effects models and extensions in ecology with R. Springer, New York.