Kootenay Lake Exploitation Study 2013

Data

The data were collected by releasing large (\(\ge\) 500 mm) rod caught rainbow trout and bull trout with acoustic tags and T-bar reward tags. An acoustic receiver array detected movements of the fish while anglers reported recaptures.

The back-calculation of length-at-age from scale samples was conducted by Greg Andrusak of Redfish Consulting Ltd. The Ministry of Forests, Lands and Natural Resource Operations (MFLNRO) provided spawner escapement estimates for kokanee in Meadow Creek spawning channel and rainbow trout at Gerrard as well as data from the Kootenay Lake Rainbow Trout (KLRT) mailout survey.

Statistical Analysis

Hierarchical Bayesian models were fitted to the rainbow and bull trout tagging, acoustic detection and angler return data for Kootenay Lake using R version 3.0.3 (Team, 2013) and JAGS 3.3.0 (Plummer, 2012) which interfaced with each other via jaggernaut 1.7 (Thorley, 2014). 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) pages 41-44.

Unless specified, the models assumed vague (low information) prior distributions (Kéry and Schaub, 2011, p. 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 (Kéry and Schaub, 2011, pp. 38-40). Model convergence was confirmed by ensuring that Rhat (Kéry and Schaub, 2011, p. 40) was less than 1.1 for each of the parameters in the model (Kéry and Schaub, 2011, p. 61). When possible model adequacy was confirmed by examination of residual plots.

The posterior distributions of the fixed (Kéry and Schaub, 2011, p. 75) parameters are summarised below 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 (Kéry and Schaub, 2011, p. 37,42).

When none of the candidate models included random effects, model selection was achieved by choosing the model with the lower Deviance Information Criterion (DIC) value. Otherwise model selection was achieved by percent relative error-based backwards stepwise variable selection (Kéry and Schaub, 2011, p. 42) from the full model until all the variables had a percent relative error less than 100. This is approximately equivalent to dropping insignificant (p > 0.05) fixed variables and uninformative random (Kéry and Schaub, 2011, pp. 77-82) variables.

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) (Kéry and Schaub, 2011, pp. 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% credible intervals (Bradford et al. 2005). Through the report bull trout data and estimates are plotted in black while rainbow trout are plotted in red. Plots were produced using the ggplot2 R package (Wickham, 2009).

Body Condition

Temporal variation in condition (body weight when accounting for body length) was estimated from the initial captures using a mass-length model (He et al. 2008).

Key assumptions of the full condition model include:

• Weight varies with fork length and day of the year (as a second order polynomial).
• Weight varies randomly with year.
• The residual variation in weight is log-normally distributed.

Competing models were fitted without year and/or day of the year.

Short-Term Hooking Mortality

Short-term (pre-release) hooking mortality was estimated from the initial captures using a two-way ANOVA (Kéry and Schaub, 2011, pp. 100-106) with an interaction and a binomial (logistic) response.

Key assumptions of the hooking mortality model include:

• Hooking mortality varies with fish size (small versus large) and hook type (single versus treble) and the interaction between the two.
• Observed hooking mortality is Bernoulli distributed.

Competing models were fitted without the fish size and/or hook type effects and their interaction.

T-bar Tag Loss

T-bar tag loss was estimated from the number of tags reported from recaught double-tagged individuals (Fabrizio et al. 1999).

Key assumptions of the tag loss model include:

• The probability of tag loss is independent between tags on a fish.
• Reported tag numbers are described by a zero-truncated binomial distribution (as fish that had lost both tags were not reportable).

Preliminary analyses indicated that there were insufficient recaptures to include days at large or fish length as explanatory variables.

Harvest Rates

Harvest rates were estimated from the number of fish reported caught and harvested in the KLRT mailout survey since 1999 using a logistic ANOVA.

Key assumptions of the full harvest model include:

• The harvest rate varies with weight class (2-5, 5-7 and >7kg).
• The number of fish harvested is described by an overdispersed binomial distribution.

Competing models were fitted without the weight effect.

Annual Growth

Annual growth was estimated from the initial capture’s scale-based back-calculated lengths using the Fabens form of the Von Bertalanffy growth model (He and Bence, 2007).

Key assumptions of the full growth model include:

• Annual length increments from age 3 and older are described by a Von Bertalanffy growth curve.
• The length-at-infinity is 1000 mm.
• Time at length 0 is 0.
• The growth coefficient (k) varies randomly with year and fish.
• Residual annual length increments are normally distributed.

Competing models were fitted without the annual and/or individual variation in k.

Natural and Fishing Mortality

Natural and fishing mortality was estimated from the acoustic detection and tag reporting data using a Cormack-Jolly-Seber state-space model (Kéry and Schaub, 2011, pp. 175-177).

Key assumptions of the mortality model include:

• One year consists of four discrete time intervals (seasons) of 3 months.
• Log-odds probability of the loss of a T-bar tag is normally distributed with a mean of -1.37 and a SD of 0.40 for bull trout and a mean of -2.75 and SD of 0.57 for rainbow trout as estimated by the T-bar Tag Loss models.
• Non-reporting of T-bar tagged angled fish uniformly ranges between 0 and 10%.
• Change in length is as described by a Von Bertalanffy growth curve with a length-at-infinity of 1000 mm and a growth coefficient (k) of 0.19.
• Probability of spawning is independent of whether a fish spawned in the previous year.
• Probability of spawning varies with length.
• Probability of capture by an angler varies randomly with year.
• Natural mortality differs between in-lake versus spawners.
• In-lake natural mortality varies with length.

Preliminary analyses indicated that length and season were not significant predictors and year was not an informative predictor of capture probability.

Model Code

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

model {

  bWeight ~ dnorm(0, 2^-2)
  bWeightLength ~ dnorm(3, 2^-2)
  bWeightDayte ~ dnorm(0, 2^-2)
  bWeightDayte2 ~ dnorm(0, 2^-2)

  sWeightYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bWeightYear[i] ~ dnorm(0, sWeightYear^-2)
  }

  sWeight ~ dunif(0, 2)

  for(i in 1:length(Length)) {

    log(eWeight[i]) <- bWeight 
        + bWeightLength * Length[i] 
        + bWeightDayte * Dayte[i] 
        + bWeightDayte2 * Dayte[i]^2 
        + bWeightYear[Year[i]]

    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }
}

model {

  bWeight ~ dnorm(0, 2^-2)
  bWeightLength ~ dnorm(3, 2^-2)

  sWeightYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bWeightYear[i] ~ dnorm(0, sWeightYear^-2)
  }

  sWeight ~ dunif(0, 2)

  for(i in 1:length(Length)) {

    log(eWeight[i]) <- bWeight 
      + bWeightLength * Length[i] 
      + bWeightYear[Year[i]]

    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }
}

model {

  bWeight ~ dnorm(0, 2^-2)
  bWeightLength ~ dnorm(3, 2^-2)
  bWeightDayte ~ dnorm(0, 2^-2)
  bWeightDayte2 ~ dnorm(0, 2^-2)

  sWeight ~ dunif(0, 2)

  for(i in 1:length(Length)) {

    log(eWeight[i]) <- bWeight 
      + bWeightLength * Length[i] 
      + bWeightDayte * Dayte[i] 
      + bWeightDayte2 * Dayte[i]^2 

    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }
}

model {

  bWeight ~ dnorm(0, 2^-2)
  bWeightLength ~ dnorm(3, 2^-2)

  sWeight ~ dunif(0, 2)

  for(i in 1:length(Length)) {

    log(eWeight[i]) <- bWeight 
      + bWeightLength * Length[i] 

    Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
  }
}

model{

  bMortality ~ dnorm(0, 5^-2)

  bMortalitySize[1] <- 0
  for (i in 2:nSize) {
    bMortalitySize[i] ~ dnorm(0, 5^-2)
  }

  bMortalityHookType[1] <- 0
  for (i in 2:nHookType) {
    bMortalityHookType[i] ~ dnorm(0, 5^-2)
  }

  for (j in 1:nHookType) {
    bMortalitySizeHookType[1, j] <- 0
  }

  for (i in 2:nSize) {
      bMortalitySizeHookType[i, 1] <- 0
      for (j in 2:nHookType) {
          bMortalitySizeHookType[i, j] ~ dnorm(0, 5^-2)
      }
  }

  for(i in 1:length(Mortality)) {

    logit(eMortality[i]) <- bMortality 
      + bMortalitySize[Size[i]] 
      + bMortalityHookType[HookType[i]]
      + bMortalitySizeHookType[Size[i], HookType[i]]

    Mortality[i] ~ dbern(eMortality[i])
  }
}

model{

  bMortality ~ dnorm(0, 5^-2)

  bMortalitySize[1] <- 0
  for (i in 2:nSize) {
    bMortalitySize[i] ~ dnorm(0, 5^-2)
  }

  bMortalityHookType[1] <- 0
  for (i in 2:nHookType) {
    bMortalityHookType[i] ~ dnorm(0, 5^-2)
  }

  for(i in 1:length(Mortality)) {

    logit(eMortality[i]) <- bMortality 
      + bMortalitySize[Size[i]] 
      + bMortalityHookType[HookType[i]]

    Mortality[i] ~ dbern(eMortality[i])
  }
}

model{

  bMortality ~ dnorm(0, 5^-2)

  bMortalitySize[1] <- 0
  for (i in 2:nSize) {
    bMortalitySize[i] ~ dnorm(0, 5^-2)
  }

  for(i in 1:length(Mortality)) {

  logit(eMortality[i]) <- bMortality 
    + bMortalitySize[Size[i]] 

  Mortality[i] ~ dbern(eMortality[i])
  }
}

model{

  bMortality ~ dnorm(0, 5^-2)

  bMortalityHookType[1] <- 0
  for (i in 2:nHookType) {
    bMortalityHookType[i] ~ dnorm(0, 5^-2)
  }

  for(i in 1:length(Mortality)) {

    logit(eMortality[i]) <- bMortality 
      + bMortalityHookType[HookType[i]]

    Mortality[i] ~ dbern(eMortality[i])
  }
}

model{

  bMortality ~ dnorm(0, 5^-2)

  for(i in 1:length(Mortality)) {

  logit(eMortality[i]) <- bMortality 

  Mortality[i] ~ dbern(eMortality[i])
  }
}

model{
  bTagLossIntercept ~ dnorm(0, 2^-2)

  for (i in 1:length(TagsRecap)) {
    logit(eTagLoss[i]) <- bTagLossIntercept

    TagsRecap[i] ~ dbin(1 - eTagLoss[i], 2) T(1, )
  }
}

model{

  bHarvest ~ dnorm(0, 2^-2)

  bHarvestWeight[1] <- 0
  for(i in 2:nWeight) {
    bHarvestWeight[i] ~ dnorm(0, 2^-2)
  }

  sDispersion ~ dunif(0, 2)
  for(i in 1:length(Harvest)) {

    eDispersion[i] ~ dnorm(0, sDispersion^-2)

    logit(eHarvest[i]) <- bHarvest
      + bHarvestWeight[Weight[i]]
      + eDispersion[i]

    Harvest[i] ~ dbin(eHarvest[i], Catch[i])
  }
}

model{

  bHarvest ~ dnorm(0, 2^-2)

  sDispersion ~ dunif(0, 2)
  for(i in 1:length(Harvest)) {

    eDispersion[i] ~ dnorm(0, sDispersion^-2)

    logit(eHarvest[i]) <- bHarvest
      + eDispersion[i]

    Harvest[i] ~ dbin(eHarvest[i], Catch[i])
  }
}

model{

  kIntercept ~ dunif (0, 1)
  Linf <- 1000

  sKYear ~ dunif(0, 1)
  for (i in 1:nYear) {
    bKYear[i] ~ dnorm(0, sKYear^-2)
  }

  sKFish ~ dunif(0, 1)
  for (i in 1:nFish) {
    bKFish[i] ~ dnorm(0, sKFish^-2)
  }

  sGrowth ~ dunif(0, 100)

  for (i in 1:length(Growth)) {

    eK[i] <- kIntercept
      + bKYear[Year[i]] 
      + bKFish[Fish[i]]

    eGrowth[i] <- (Linf - Length[i]) * (1 - exp(-eK[i]))
    Growth[i] ~ dnorm(eGrowth[i], sGrowth^-2)
  }
}

model{

  kIntercept ~ dunif (0, 1)
  Linf <- 1000

  sKYear ~ dunif(0, 1)
  for (i in 1:nYear) {
    bKYear[i] ~ dnorm(0, sKYear^-2)
  }

  sGrowth ~ dunif(0, 100)

  for (i in 1:length(Growth)) {

  eK[i] <- kIntercept
    + bKYear[Year[i]] 

  eGrowth[i] <- (Linf - Length[i]) * (1 - exp(-eK[i]))
  Growth[i] ~ dnorm(eGrowth[i], sGrowth^-2)
  }
}

model {

  kIntercept ~ dunif (0, 1)
  Linf <- 1000

  sKFish ~ dunif(0, 1)
  for (i in 1:nFish) {
    bKFish[i] ~ dnorm(0, sKFish^-2)
  }

  sGrowth ~ dunif(0, 100)

  for (i in 1:length(Growth)) {

    eK[i] <- kIntercept
    + bKFish[Fish[i]]

    eGrowth[i] <- (Linf - Length[i]) * (1 - exp(-eK[i]))
    Growth[i] ~ dnorm(eGrowth[i], sGrowth^-2)
  }
}

model {

  kIntercept ~ dunif (0, 1)
  Linf <- 1000

  sGrowth ~ dunif(0, 100)

  for (i in 1:length(Growth)) {

    eK[i] <- kIntercept

    eGrowth[i] <- (Linf - Length[i]) * (1 - exp(-eK[i]))
    Growth[i] ~ dnorm(eGrowth[i], sGrowth^-2)
  }
}

model {

  bHandlingM ~ dnorm(0, 2^-2)
  bAngled ~ dnorm(0, 2^-2)
  bHarvest ~ dnorm(0, 2^-2)
  bNaturalM ~ dnorm(0, 2^-2)
  bSpawned ~ dnorm(0, 2^-2)

  bTagLoss ~ dnorm(TagLoss[1], TagLoss[2]^-2)
  logit(pTagLoss[1]) <- bTagLoss
  for (i in 2:nTags) {
    pTagLoss[i] <- pTagLoss[1]^2
  }

  pNonReporting ~ dunif(NonReporting[1], NonReporting[2])

  bNaturalMSpawning ~ dnorm(0, 2^-2)

  bSpawnedLength ~ dnorm(0, 2^-2)
  bNaturalMLength ~ dnorm(0, 2^-2)

  for (i in 1:nFish) {

    logit(eHandlingM[i, Monitored[i]]) <- bHandlingM
    logit(eAngled[i, Monitored[i]]) <- bAngled 
    logit(eHarvest[i, Monitored[i]]) <- bHarvest
    logit(eSpawned[i, Monitored[i]]) <- bSpawned 
                                      + bSpawnedLength * Length[i, Monitored[i]]

    dAliveM[i] ~ dbern((1-eHandlingM[i, Monitored[i]]))

    Spawned[i, Monitored[i]] ~ dbern(dAliveM[i] * eSpawned[i, Monitored[i]]
                                     * equals(Season[Monitored[i]], SpawningSeason))

    logit(eNaturalM[i, Monitored[i]]) <- 
                        (bNaturalM + bNaturalMLength * Length[i, Monitored[i]]) 
                            * (1-Spawned[i, Monitored[i]])
                        + bNaturalMSpawning * Spawned[i, Monitored[i]]

    dAngled[i, Monitored[i]] ~ dbern(dAliveM[i] * eAngled[i, Monitored[i]])

    Reported[i, Monitored[i]] ~ dbern(
      dAngled[i, Monitored[i]] * (1 - pTagLoss[Tags[i]]) * (1 - pNonReporting)
    )

    Harvested[i, Monitored[i]] ~ dbern(
      dAngled[i, Monitored[i]] * eHarvest[i, Monitored[i]]
    )

    Alive[i, Monitored[i]] ~ dbern(
      dAliveM[i] * (1-eNaturalM[i, Monitored[i]]) 
        * (1 - (dAngled[i, Monitored[i]] * Harvested[i, Monitored[i]]
        + dAngled[i, Monitored[i]] * (1-Harvested[i, Monitored[i]]) 
            * eHandlingM[i, Monitored[i]]))
    )

    for (j in (Monitored[i]+1):nInterval) {

      logit(eHandlingM[i, j]) <- bHandlingM
      logit(eAngled[i, j]) <- bAngled

      logit(eHarvest[i, j]) <- bHarvest
      logit(eSpawned[i, j]) <- bSpawned
                                  + bSpawnedLength * Length[i, j]

      Spawned[i, j] ~ dbern(
        Alive[i, j-1] * eSpawned[i, j] * equals(Season[j],SpawningSeason)
      )

      logit(eNaturalM[i, j]) <- (bNaturalM + bNaturalMLength * Length[i, j]) 
                                  * (1-Spawned[i, j])
                                + bNaturalMSpawning * Spawned[i, j]

      dAngled[i, j] ~ dbern(Alive[i, j-1] * eAngled[i, j])

      Reported[i, j] ~ dbern(
          dAngled[i, j] * (1 - pTagLoss[Tags[i]]) * (1 - pNonReporting)
      )

    Harvested[i, j] ~ dbern(
      dAngled[i, j] * eHarvest[i, j]
    )
      Alive[i,j] ~ dbern(Alive[i, j-1] * (1-eNaturalM[i,j]) 
                         * (1 - (dAngled[i, j] * Harvested[i, j]
                                 + dAngled[i, j] * (1-Harvested[i, j]) 
                                 * eHandlingM[i, j])))
    }
  }
}

Model Parameters

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

Figures

Escapement

Fishery

Age - Rainbow Trout

Detections

Body Condition

Short-Term Hooking Mortality

T-Bar Tag Loss

Harvest Rates

Growth - Rainbow Trout

Natural And Fishing Mortality