# 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: http://www.poissonconsulting.ca/f/1007329429.

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

## Background

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:

- Estimate the observer efficiency when conducting spring snorkel surveys for age-1 Rainbow Trout.
- Estimate the spring abundance of age-1 Rainbow Trout by year.
- 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.

## Methods

### 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:

- The observer probability varies with study design (Index versus GPS).
- There is no tag loss.
- There is no emigration of marked fish.
- The number of marked fish that are resighted is described by a binomial distribution.

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

#### Abundance

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:

- The lineal fish density varies with year, useable width and river kilometer as a polynomial, and randomly with site.
- The observer efficiency at marked sites was as estimated by the observer efficiency model (which varied by study design).
- The observer efficiency also varies with visit type (standard count versus presence of marked fish) within study design, and randomly with swimmer.
- The expected count at a site is the expected lineal density multiplied by the site length, the observer efficiency and the proportion of the site surveyed.
- The residual variation in the actual count is gamma-Poisson distributed.

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

#### Stock-Recruitment

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 prior probability \(\alpha\) is normally distributed with a mean of 500 and a SD of 250; this mean is based on an average of 8,000 eggs per female spawner, a 50:50 sex ratio, 50% egg survival, 50% post-emergence fall survival and 50% overwintering survival.
- The residual variation in the number of recruits is log-normally distributed.

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

#### Escapement

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 `ii` ^{th} study design on `logit(eEfficiency)` |

`eEfficiency[ii]` |
Expected capture efficiency on `ii` ^{th} visit |

`Marked[ii]` |
Number of marked fish prior to `ii` ^{th} visit |

`Resighted[ii]` |
Number of marked fish resighted on `ii` ^{th} visit |

`StudyDesign[ii]` |
Study design of `ii` ^{th} visit |

##### Capture Efficiency - Model1

```
model{
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])
}
}
```

#### Abundance

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 `ii` ^{th} site on `log(eDensity)` |

`bDensityWidth` |
Exponent for effect of width on `eDensity` |

`bDensityYear[ii]` |
Effect of `ii` ^{th} year on `log(eDensity)` |

`bEfficiencySwimmer[ii]` |
Effect of `ii` ^{th} swimmer on `logit(eEfficiency)` |

`bEfficiencyVisitStudy[ii, jj]` |
Effect of `ii` ^{th} visit type within `jj` ^{th} study design on `logit(eEfficiency)` |

`Count[ii]` |
Total number of fish observed on `ii` ^{th} visit |

`eAbundance[ii]` |
Expected abundance of fish at site of `ii` ^{th} visit |

`eCount[ii]` |
Expected total number of fish at site of `ii` ^{th} visit |

`eDensity[ii]` |
Expected lineal density of fish at site of `ii` ^{th} visit |

`eDispersion` |
Overdispersion of `Count` |

`eEfficiency[ii]` |
Expected observer efficiency on `ii` ^{th} visit |

`logit(bEfficiencyStudy[ii])` |
Effect of `ii` ^{th} 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 `ii` ^{th} 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 `ii` ^{th} visit |

`SiteLength[ii]` |
Length of site of `ii` ^{th} visit |

`StudyDesign[ii]` |
Study design of `ii` ^{th} visit |

`SurveyProportion[ii]` |
Proportion of site surveyed on `ii` ^{th} visit |

`Swimmer[ii]` |
Swimmer of `ii` ^{th} visit |

`VisitType[ii]` |
Visit type of `ii` ^{th} visit |

`Width[ii]` |
Site width of `ii` ^{th} visit |

`Year[ii]` |
Year of `ii` ^{th} visit |

##### Abundance - Model1

```
model{
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], 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])
}
}
```

#### Stock-Recruitment

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

```
model{
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)
}
}
```

#### Escapement

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

```
model{
bRinf ~ dnorm(Rinf.mean, Rinf.sd^-2) T(Rinf.lower, Rinf.upper)
bAlpha ~ dnorm(alpha.mean, alpha.sd^-2) T(alpha.lower, alpha.upper)
bS0 <- bRinf/bAlpha
}
```

## Results

### 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 |

#### Stock-Recruitment

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 |

#### Escapement

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 |

### Figures

#### Capture Efficiency - Age-1

#### Abundance - Age-1

#### Stock-Recruitment

#### Escapement

## Acknowledgements

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

- Habitat Conservation Trust Foundation (HCTF) and the anglers, hunters, trappers and guides who contribute to the Trust
- Fish and Wildlife Compensation Program (FWCP) and its program partners BC Hydro, the Province of BC and Fisheries and Oceans Canada
- Ministry of Forests, Lands and Natural Resource Operations (MFLNRO)
- Redfish Consulting
- Greg Andrusak

- Gary Pavan

## References

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. http://a100.gov.bc.ca/pub/acat/public/viewReport.do?reportId=36011.

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. https://doi.org/10.1139/f05-179.

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. https://doi.org/10.1111/j.1439-0426.2008.01187.x.

Kery, Marc, and Michael Schaub. 2011. *Bayesian Population Analysis Using WinBUGS : A Hierarchical Perspective*. Boston: Academic Press. http://www.vogelwarte.ch/bpa.html.

Plummer, Martyn. 2012. “JAGS Version 3.3.0 User Manual.” http://sourceforge.net/projects/mcmc-jags/files/Manuals/3.x/.

Team, R Core. 2013. “R: A Language and Environment for Statistical Computing.” Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org.

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