Lardeau-Lower Duncan River Age-1 Rainbow Trout Abundance and Stock-Recruitment Analysis

Suggested Citation: Thorley, J.L. and Hogan, P.M. (2015) Lardeau-Lower Duncan River Age-1 Rainbow Trout Abundance and Stock-Recruitment Analysis 2015. A Poisson Consulting Analysis Report. URL:

The interpretive report from 2013 (Andrusak 2013) is available online.


Juvenile Gerrard Rainbow Trout from Kootenay Lake rear in the Lardeau and Lower Duncan rivers. Between 2006 and 2014 spring snorkel surveys have been undertaken to estimate the abundance of the age-1 fish. Between 2006 and 2010 the surveys were conducted at fixed index site. From 2010 fish observations were mapped to sites on the river by georeferencing all counts which allowed the swimmers to cover more of the river.

The primary aims of the current analyses are to:

  1. Estimate the observer efficiency when conducting spring snorkel surveys for age-1 Rainbow Trout.
  2. Estimate the spring abundance of age-1 Rainbow Trout by year.
  3. Estimate the stock-recruitment relationship between the number of Rainbow Trout spawners in a given year and the number of age-1 recruits the following spring.


Data Preparation

The snorkel count data was provided by Redfish Consulting Ltd. The Ministry of Forests, Lands and Natural Resource Operations (MFLNRO) provided the AUC-based spawner escapement estimates for Gerrard.

Length-Bias Correction

Biases in length estimation were estimated manually; for each observer for each unique day of observations, a fork-length cut-off between age-1 and age-2+ Rainbow Trout was estimated, which was then used to classify each fish by lifestage.

Statistical Analysis

Hierarchical Bayesian models were fitted to the data using R version 3.2.2 (Team 2013) and JAGS 3.4.0 (Plummer 2012) which interfaced with each other via jaggernaut 2.3.1 (Thorley 2013). For additional information on hierarchical Bayesian modelling in the BUGS language, of which JAGS uses a dialect, the reader is referred to (Kery and Schaub 2011, 41–44).

Unless indicated otherwise, the models use prior distributions that are vague in the sense that they did not affect the posterior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from a minimum of 1,000 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of three chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that Rhat (Kery and Schaub 2011, 40) was less than 1.1 for each of the parameters in the model (Kery and Schaub 2011, 61).

The posterior distributions of the fixed (Kery and Schaub 2011, 75) parameters are summarised in terms of a point estimate (mean), lower and upper 95% credible limits (2.5th and 97.5th percentiles), the standard deviation (SD), percent relative error (half the 95% credible interval as a percent of the point estimate) and significance (Kery and Schaub 2011, 37, 42).

Variable selection was achieved by dropping insignificant (Kery and Schaub 2011, 37, 42) fixed (Kery and Schaub 2011, 77–82) variables and uninformative random variables. A fixed variables was considered to be insignificant if its significance was \(\geq\) 0.05 while a random variable was considered to be uninformative if its percent relative error was \(\geq\) 80%. The Deviance Information Criterion (DIC) was not used because it is of questionable validity when applied to hierarchical models (Kery and Schaub 2011, 469).

The results are displayed graphically by plotting the modeled relationships between particular variables and the response with 95% credible intervals (CRIs) with the remaining variables held constant. In general, continuous and discrete fixed variables are held constant at their mean and first level values respectively while random variables are held constant at their typical values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). Where informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% CRIs (Bradford, Korman, and Higgins 2005).

Observer Efficiency

The observer efficiency for age-1 Rainbow Trout was estimated using a mark-resight binomial model (Kery and Schaub 2011, 134–36, 384–88).

Key assumptions of the observer efficiency model include:

Useable width and swimmer were not found to be significant or informative predictors of observer efficiency, respectively.


The abundance was estimated from the bias-corrected observer count data using an overdispersed Poisson model (Kery and Schaub 2011, 55–56). The annual abundance estimates represent the total number of fish in the study area.

Key assumptions of the abundance model include:

Curvature and channel class were not found to be significant predictors of density.


The relationship between the number of Rainbow Trout spawners in a given year (\(S\)) and the number of age-1 Rainbow Trout the following spring (\(R\)) was estimated using the smooth hockey stick stock-recuitment model (Froese 2008):

\[ R = R_\infty \cdot \left(1 - \exp{\left\{-\frac{\alpha}{R_\infty} \cdot S\right\}} \right) \quad,\]

where \(R_\infty\) is the carrying capacity of recruits, and \(\alpha\) is the number of recruits per spawner at low density.

Key assumptions of the stock-recruitment model include:

The Q90 discharge in July and the Q10 discharge from January to February were not found to be significant predictors of \(R_\infty\).


The number of spawners required to reach a given percentage the recruit carrying capacity (\(\epsilon = R/R_\infty\)) was estimated from the hockey stick stock-recuitment model using the equation:

\[ S = \frac{R_\infty}{\alpha} \cdot \ln{\left(1 - \epsilon\right)^{-1}} \quad. \]

Model Code

The JAGS model code, which uses a series of naming conventions, is presented below.

Capture Efficiency

Variable/Parameter Description
bEfficiency logit(eEfficiency) intercept
bEfficiencyStudyDesign[ii] Effect of iith study design on logit(eEfficiency)
eEfficiency[ii] Expected capture efficiency on iith visit
Marked[ii] Number of marked fish prior to iith visit
Resighted[ii] Number of marked fish resighted on iith visit
StudyDesign[ii] Study design of iith visit
Capture Efficiency - Model1

  bEfficiency ~ dnorm(0, 5^-2)

  bEfficiencyStudyDesign[1] <- 0
  for (ii in 2:nStudyDesign) {
    bEfficiencyStudyDesign[ii] ~ dnorm(0, 2^-2)

  for(ii in 1:length(Marked)){
    logit(eEfficiency[ii]) <- bEfficiency
                            + bEfficiencyStudyDesign[StudyDesign[ii]]

    Resighted[ii] ~ dbin(eEfficiency[ii], Marked[ii])


Variable/Parameter Description
bDensityMarking Effect of Marking on log(eDensity)
bDensityRkmX Polynomial coefficients of effect of river kilometer on log(eDensity)
bDensitySite[ii] Effect of iith site on log(eDensity)
bDensityWidth Exponent for effect of width on eDensity
bDensityYear[ii] Effect of iith year on log(eDensity)
bEfficiencySwimmer[ii] Effect of iith swimmer on logit(eEfficiency)
bEfficiencyVisitStudy[ii, jj] Effect of iith visit type within jjth study design on logit(eEfficiency)
Count[ii] Total number of fish observed on iith visit
eAbundance[ii] Expected abundance of fish at site of iith visit
eCount[ii] Expected total number of fish at site of iith visit
eDensity[ii] Expected lineal density of fish at site of iith visit
eDispersion Overdispersion of Count
eEfficiency[ii] Expected observer efficiency on iith visit
logit(bEfficiencyStudy[ii]) Effect of iith study design on logit(eEfficiency)
Marking Whether a site has been chosen as a marking site under the different study designs
Rkm[ii] River kilometer of iith visit
sDensitySite SD of effect of site on log(eDensity)
sDispersion SD of overdispersion term
sEfficiencySwimmer SD of effect of swimmer on logit(eEfficiency)
Site[ii] Site of iith visit
SiteLength[ii] Length of site of iith visit
StudyDesign[ii] Study design of iith visit
SurveyProportion[ii] Proportion of site surveyed on iith visit
Swimmer[ii] Swimmer of iith visit
VisitType[ii] Visit type of iith visit
Width[ii] Site width of iith visit
Year[ii] Year of iith visit
Abundance - Model1

  for(yr in 1:nYear){
    bDensityYear[yr] ~ dnorm(0, 5^-2)

  bDensityWidth ~ dnorm(0, 2^-2)

  sDensitySite ~ dunif(0, 5)
  for(st in 1:nSite){
    bDensitySite[st] ~ dnorm(-sDensitySite^2 / 2, sDensitySite^-2)

  bDensityMarking[1] <- 0
  for(mk in 2:nMarking) {
    bDensityMarking[mk] ~ dnorm(0, 5^-2)

  sEfficiencySwimmer ~ dunif(0, 5)
  for(sw in 1:nSwimmer){
    bEfficiencySwimmer[sw] ~ dnorm(0, sEfficiencySwimmer^-2)

  bDensityRkm1 ~ dnorm(0, 2^-2)
  bDensityRkm2 ~ dnorm(0, 2^-2)
  bDensityRkm3 ~ dnorm(0, 2^-2)
  bDensityRkm4 ~ dnorm(0, 2^-2)

  for(sd in 1:nStudyDesign){
    bEfficiencyStudy[sd] ~ dnorm(Efficiency[sd],[sd]^-2) T(Efficiency.lower[sd], Efficiency.upper[sd])

  for(sd in 1:nStudyDesign){
    bEfficiencyVisitStudy[1, sd] <- 0
    for(vt in 2:nVisitType){
      bEfficiencyVisitStudy[vt, sd] ~ dnorm(0, 2^-2)

  sDispersion ~ dunif(0, 5)

  for(ii in 1:length(Count)){
    log(eDensity[ii]) <- bDensityYear[Year[ii]]
                       + bDensityWidth * log(Width[ii])
                       + bDensitySite[Site[ii]]
                       + bDensityRkm1 * Rkm[ii]
                       + bDensityRkm2 * Rkm[ii]^2
                       + bDensityRkm3 * Rkm[ii]^3
                       + bDensityRkm4 * Rkm[ii]^4
                       + bDensityMarking[Marking[ii]]

    eAbundance[ii] <- eDensity[ii] * SiteLength[ii]

    logit(eEfficiency[ii]) <- logit(bEfficiencyStudy[StudyDesign[ii]])
                            + bEfficiencyVisitStudy[VisitType[ii], StudyDesign[ii]]
                            + bEfficiencySwimmer[Swimmer[ii]]

    eCount[ii] <- eAbundance[ii] * eEfficiency[ii] * SurveyProportion[ii]

    eDispersion[ii] ~ dgamma(1/sDispersion^2, 1/sDispersion^2)

    Count[ii] ~ dpois(eCount[ii] * eDispersion[ii])


Variable/Parameter Description
alpha Recruits per spawner at low density
eRecruits Expected number of recruits
Recruits Number of recruits
Rinf Carrying capacity of recruits
Spawners Number of spawners
sRecruits SD of Recruits
Stock-Recruitment - Model1

  alpha ~ dnorm(8000 * 0.5^4, 250^-2) T(0,)
  Rinf ~ dunif(5 * 10^4, 2.5 * 10^5)

  sRecruits ~ dunif(0, 5)
  for (ii in 1:length(Spawners)) {
    eRecruits[ii] <- Rinf * (1 - exp(-alpha / Rinf *  Spawners[ii]))
    Recruits[ii] ~ dlnorm(log(eRecruits[ii]) - sRecruits^2 / 2, sRecruits^-2)


Variable/Parameter Description
bAlpha Expected number of recruits per spawner at low density
bRinf Expected carrying capacity of recruits
bS0 Predicted spawner coefficient
Escapement - Model1

  bRinf ~ dnorm(Rinf.mean,^-2) T(Rinf.lower, Rinf.upper)
  bAlpha ~ dnorm(alpha.mean,^-2) T(alpha.lower, alpha.upper)

  bS0 <- bRinf/bAlpha


Model Parameters

The posterior distributions for the fixed (Kery and Schaub 2011 p. 75) parameters in each model are summarised below.

Capture Efficiency - Age-1

Parameter Estimate Lower Upper SD Error Significance
bEfficiency 0.2220 -0.1006 0.5289 0.1579 140 0.1639
bEfficiencyStudyDesign[2] -0.9689 -1.3826 -0.5597 0.2118 42 0.0010
Convergence Iterations
1 5000

Abundance - Age-1

Parameter Estimate Lower Upper SD Error Significance
bDensityMarking[2] 0.3541 0.0283 0.6912 0.1704 94 0.0323
bDensityMarking[3] 0.3187 0.0884 0.5435 0.1195 71 0.0133
bDensityRkm1 -0.2544 -0.4290 -0.0770 0.0891 69 0.0010
bDensityRkm2 0.1748 -0.0461 0.4127 0.1198 130 0.1425
bDensityRkm3 -0.0555 -0.1358 0.0259 0.0398 150 0.1577
bDensityRkm4 -0.1479 -0.2271 -0.0713 0.0403 53 0.0010
bDensityWidth 0.2481 0.1386 0.3542 0.0532 43 0.0010
bDensityYear[1] -0.3910 -1.0980 0.3630 0.3600 190 0.2755
bDensityYear[2] -1.3740 -2.1770 -0.5830 0.4200 58 0.0019
bDensityYear[3] -1.1600 -1.9250 -0.4070 0.3840 65 0.0038
bDensityYear[4] -1.4410 -2.2550 -0.5920 0.4180 58 0.0010
bDensityYear[5] -1.3990 -2.2700 -0.6090 0.4240 59 0.0038
bDensityYear[6] -1.1366 -1.5930 -0.6249 0.2482 43 0.0010
bDensityYear[7] -1.1724 -1.6130 -0.6672 0.2491 40 0.0010
bDensityYear[8] -1.5070 -1.9234 -1.0317 0.2291 30 0.0010
bDensityYear[9] -1.9347 -2.3619 -1.4749 0.2316 23 0.0010
bEfficiencyStudy[1] 0.5542 0.4874 0.6174 0.0347 12 0.0010
bEfficiencyStudy[2] 0.3254 0.2724 0.3726 0.0263 15 0.0010
bEfficiencyVisitStudy[2,1] -1.2220 -2.1670 -0.1350 0.4970 83 0.0247
bEfficiencyVisitStudy[2,2] -0.3607 -0.7559 0.0176 0.1995 110 0.0608
sDensitySite 0.5654 0.4567 0.6619 0.0519 18 0.0010
sDispersion 0.8567 0.7896 0.9248 0.0339 8 0.0010
sEfficiencySwimmer 0.3956 0.2217 0.7127 0.1275 62 0.0010
Convergence Iterations
1.06 40000


Parameter Estimate Lower Upper SD Error Significance
alpha 5.797e+02 2.478e+02 1016.800 2.011e+02 66 0.001
Rinf 9.880e+04 6.930e+04 171600.000 2.280e+04 52 0.001
sRecruits 5.125e-01 2.928e-01 0.891 1.583e-01 58 0.001
Convergence Iterations
1.04 1e+06


Parameter Estimate Lower Upper SD Error Significance
bAlpha 592.1 286.1 935.9 171.2 55 7e-04
bRinf 103150.0 72420.0 144450.0 19340.0 35 7e-04
bS0 191.7 93.4 379.1 74.8 74 7e-04
Convergence Iterations
1 1000


Capture Efficiency - Age-1

Figure 1. Predicted capture efficiency for Age-1 Rainbow Trout by study design (with 95% CRIs).

Abundance - Age-1

Figure 2. Predicted lineal density of Age-1 Rainbow Trout by site marking type (with 95% CRIs).
Figure 3. Predicted observer efficiency for Age-1 Rainbow Trout by visit type and study design (with 95% CRIs).
Figure 4. Predicted lineal density of Age-1 Rainbow Trout in 2014 by useable width (with 95% CRIs).
Figure 5. Predicted lineal density of Age-1 Rainbow Trout by river kilometre (with 95% CRIs).
Figure 6. Predicted lineal density of Age-1 Rainbow Trout by year (with 95% CRIs).
Figure 7. Predicted abundance of Age-1 Rainbow Trout by year (with 95% CRIs).
Figure 8. Predicted abundance of Age-1 Rainbow Trout by river and year (with 95% CRIs).


Figure 9. Predicted Rainbow Trout stock-recruitment relationship between AUC spawners and age-1 recruits (with 95% CRIs).


Figure 10. Predicted number of AUC spawners required to reach a given percentage of age-1 recruit carrying capacity (with 95% CRIs).


The organisations and individuals whose contributions have made this analysis report possible include:


Andrusak, G.F. 2013. “Determination of Gerrard Rainbow Trout Productivity and Capacity in Defining Management Reference Points: Progress Report.” A Redfish Consulting Ltd. Report. Nelson, B.C.: Fish; Wildlife Compensation Program - Columbia Basin, Habitat Conservation Trust Foundation; the Ministry of Forests, Lands; Natural Resource Operations.

Bradford, Michael J, Josh Korman, and Paul S Higgins. 2005. “Using Confidence Intervals to Estimate the Response of Salmon Populations (Oncorhynchus Spp.) to Experimental Habitat Alterations.” Canadian Journal of Fisheries and Aquatic Sciences 62 (12): 2716–26.

Froese, R. 2008. “The Continuous Smooth Hockey Stick: A Newly Proposed Spawner-Recruitment Model: The New Stock-Recruitment Model.” Journal of Applied Ichthyology 24 (6): 703–4.

Kery, Marc, and Michael Schaub. 2011. Bayesian Population Analysis Using WinBUGS : A Hierarchical Perspective. Boston: Academic Press.

Plummer, Martyn. 2012. “JAGS Version 3.3.0 User Manual.”

Team, R Core. 2013. “R: A Language and Environment for Statistical Computing.” Vienna, Austria: R Foundation for Statistical Computing.

Thorley, J. L. 2013. “Jaggernaut: An R Package to Facilitate Bayesian Analyses Using JAGS (Just Another Gibbs Sampler).” Nelson, B.C.: Poisson Consulting Ltd.

blog comments powered by Disqus