Cutthroat Trout

Cowichan Lake Cuththroat Trout Exploitation Analysis 2023

Dalgarno, S. & Thorley, J.L. (2024) Cowichan Lake Cuththroat Trout Exploitation Analysis 2023. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1638637752.

Background

Cowichan Lake supports an important Cutthroat Trout fishery. To provide information on the natural and fishing mortality, Cutthroat Trout were caught by angling and tagged with acoustic transmitters and/or a $100 reward tag. Some acoustic transmitters were equipped with depth and temperature sensors.

A fixed receiver array was deployed in the lake to monitor fish movements and location year-round. Additional fixed receivers were deployed during the spawning season at known spawning tributaries. A lake census was completed 4 times per year using mobile receivers.

The objectives of this analysis are to:

Estimate fishing and natural mortality.
Estimate annual variation in fishing and natural mortality.
Estimate seasonal variation in fishing and natural mortality.
Estimate effect of handling on natural mortality.
Estimate variation in fishing and natural mortality by fork length.
Estimate spawning probability.
Estimate variation in spawning probability by fork length.
Estimate spawning mortality.

Data Preparation

The data were provided to Poisson Consulting Ltd. by Aswea Porter following extensive data screening.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The estimates were produced using JAGS (Plummer 2003). For additional information on Bayesian estimation the reader is referred to McElreath (2020).

Unless stated otherwise, the Bayesian analyses used weakly informative normal and half-normal prior distributions (Gelman, Simpson, and Betancourt 2017). The posterior distributions were estimated from 1500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of 3 chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that the potential scale reduction factor $\hat{R} \leq 1.05$ (Kery and Schaub 2011, 40) and the effective sample size (Brooks et al. 2011) $\textrm{ESS} \geq 150$ for each of the monitored parameters (Kery and Schaub 2011, 61).

The parameters are summarised in terms of the point estimate, lower and upper 95% credible limits (CLs) and the surprisal s-value (Greenland 2019). The estimate is the median (50th percentile) of the MCMC samples while the 95% CLs are the 2.5th and 97.5th percentiles. The s-value indicates how surprising it would be to discover that the true value of the parameter is in the opposite direction to the estimate (Greenland 2019). An s-value of $>$ 4.3 bits, which is equivalent to a significant p-value $<$ 0.05 (Kery and Schaub 2011; Greenland and Poole 2013), indicates that the surprise would be equivalent to throwing at least 4.3 heads in a row.

The condition that parameters describing the effects of secondary (nuisance) explanatory variable(s) have significant p-values was used as a model selection heuristic (Kery and Schaub 2011). Based on a similar argument, the condition that random effects have a standard deviation with a lower 95% CL $>$ 5% of the estimate was used as an additional model selection heuristic. Primary explanatory variables are evaluated based on their estimated effect sizes with 95% CLs (Bradford, Korman, and Higgins 2005)

The results are displayed graphically by plotting the modeled relationships between particular variables and the response(s) 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).

The analyses were implemented using R version 4.3.3 (R Core Team 2021) and the mbr family of packages.

Model Descriptions

Survival

The natural mortality and recapture probability were estimated using a Bayesian individual state-space survival model (Thorley and Andrusak 2017) with monthly intervals. The survival model incorporated natural, handling and spawning mortality, acoustic detections and movements from a fixed receiver array and quarterly lake census, vertical movements from acoustic tags equipped with depth sensors, FLOY tag loss, recapture, reporting and release.

The expected fork length at monthly period $t$ was estimated from length at capture ($L_{i,fi}$) plus the length increment expected under a Von Bertalanffy Growth Curve (Walters and Martell 2004) \[L_{i,t} = L_{i,fi} + (L_{∞} − L_{i,fi} )(1 − e^{(−\frac{1}{12} \cdot k(t −fi))})\] where $fi$ is the period at capture and the $\frac{1}{12}$ exponent adjusts for the fact that there are 12 periods per year.Linf was fixed to a biologically plausible value of 775 mm based on prior information (personal communication, R. Ptolemy); k of 0.23 was estimated from 88 aged Cutthroat Trout in Cowichan Lake (personal communication, Brendan Andersen).

The probability of spawning at estimated fork length $L$ is determined by

\[S = \frac{L^{S_p}}{L_s^{S_p} + L^{S_p}} \cdot {E_s}\]

where $L_s$ is the length at which 50% of fish are mature, $E_s$ is the probability of a mature fish spawning and $S_p$ is the maturity (as a function of length) power.

Additional key assumptions of the survival model include all but two of the assumptions in Thorley and Andrusak (2017) – the exceptions being 3 (no loss of FLOY tags) and 11 (anglers report all FLOY tags). The survival model also assumes that:

The monthly interval probability of FLOY tag loss is 0.001.
All recaptured fish that are released have their FLOY tags removed.
Reporting of a recaptured fish with a FLOY tag is likely greater than 90%.
The effect of handling mortality is restricted to the first two monthly periods after initial capture or recapture.
The effect of spawning mortality is restricted to the spawning period.
Acoustic tags are functional over their estimated battery lifespan and are not shed.
Recapture probability varies by month and size class.
Natural mortality varies by month, study year and size class.
The probability of detection depends on whether a fish is alive or has died near or far from a receiver and whether a census has occurred.

Study years span from October 01 to September 30 of the following year (i.e., study year 1 is 2019-10-01 to 2020-09-30). Since data collection started in October, use of study year instead of calendar year ensured that all months were represented in each year. Observations after study year 4 (i.e., in October 2023) were removed from this year’s analysis.

Size classes are defined as small ($\le$ 350 mm), medium (351 $-$ 500 mm) and large (> 500 mm). Size class membership can change in each period depending on the fork length predicted by the growth model.

Preliminary analyses indicated lack of support for variation in recapture probability and spawning probability by study year.

We assumed that the monthly interval probability of FLOY tag loss was 0.001 (0.012 annual probability) based on tag loss estimates from double-tagged Lake Trout in Horse Lake.

Preliminary attempts to estimate tag loss in the model using an uninformative prior confirmed that this assumption was consistent with the available information, although we had to fix tag loss in the model to avoid poor convergence.

We assumed that fish detected in monitored spawning tributaries represents all spawning tagged fish in the lake. The ‘spawning window’ was defined as February to May given that there was consistent receiver coverage in spawning tributaries across years. Receivers were also placed in spawning tributaries in January in study years 3 and 4. Four fish that were detected in a spawning tributary only in January (all in study year 4) were not assigned a status of ‘spawning’. Two study years had partial coverage in February (starting Feb 10-11 and Feb 17-19). These spawning detections were included because we assumed that most spawning events would occur toward the end of the month.

Model Templates

Survival

.model{
  bMortality ~ dnorm(-3, 3^-2)
  bMortalityHandling ~ dnorm(0, 2^-2)
  bMonthsHandling <- 2
  bMortalitySpawning ~ dnorm(0, 2^-2)
  bLs ~  dunif(100, 1000)
  bSp ~ dexp(1/20)
  bEs ~ dbeta(1,1)

  bDetectedAlive ~ dbeta(1, 1)
  bDieNear ~ dbeta(1, 1)
  bDetectedDeadNear ~ dbeta(10, 1)
  bDetectedDeadFar ~ dbeta(1, 10)
  bDetectedCensus ~ dbeta(1, 1)
  bMovedH ~ dbeta(1, 1)
  bMovedV ~ dbeta(1, 1)
  bRecaptured ~ dnorm(0, 2^-2)
  bReported ~ dbeta(1, 1)
  bReleased ~ dbeta(1, 1)
  bTagLoss <- 0.001

  bMortalityAnnual[1] <- 0
  for (i in 2:nAnnual) {
    bMortalityAnnual[i] ~ dnorm(0, 2^-2)
  }
  bRecapturedSize[1] <- 0
  for (i in 2:nSize) {
    bRecapturedSize[i] ~ dnorm(0, 2^-2)
  }
  bMortalitySize[1] <- 0
  for (i in 2:nSize) {
    bMortalitySize[i] ~ dnorm(0, 2^-2)
  }
  sRecapturedMonth ~ dexp(1)
  for (i in 1:nMonth) {
    bRecapturedMonth[i] ~ dnorm(0, sRecapturedMonth^-2)
  }
  sMortalityMonth ~ dexp(1)
  for (i in 1:nMonth) {
    bMortalityMonth[i] ~ dnorm(0, sMortalityMonth^-2)
  }
  for (i in 1:nCapture){
    eMortalityHandling[i, PeriodCapture[i]] <- ilogit(bMortalityHandling)
    eMonthsSinceCapture[i,PeriodCapture[i]] <- 1-Recaptured[i,PeriodCapture[i]]
    eMortalitySpawning[i, PeriodCapture[i]] <- ilogit(bMortalitySpawning) * SpawningPeriod[PeriodCapture[i]]
    
    eDetectedAlive[i] <- bDetectedAlive 
    eDieNear[i] ~ dbern(bDieNear)
    eDetectedDeadNear[i] <- bDetectedDeadNear 
    eDetectedDeadFar[i] <- bDetectedDeadFar 
    eDetectedDead[i] <- eDieNear[i] * eDetectedDeadNear[i] + 
      (1 - eDieNear[i]) * eDetectedDeadFar[i]
    
    logit(eMortality[i,PeriodCapture[i]]) <- bMortality + bMortalityAnnual[Annual[PeriodCapture[i]]] + bMortalitySize[Size[i, PeriodCapture[i]]] + bMortalityMonth[Month[PeriodCapture[i]]]

    eMovedH[i,PeriodCapture[i]] <- bMovedH
    eMovedV[i,PeriodCapture[i]] <- bMovedV
    eDetected[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * eDetectedAlive[i] +
      (1-Alive[i,PeriodCapture[i]]) * eDetectedDead[i]
    eDetectedCensus[i,PeriodCapture[i]] <- bDetectedCensus
    logit(eRecaptured[i,PeriodCapture[i]]) <- bRecaptured + bRecapturedMonth[Month[PeriodCapture[i]]] + bRecapturedSize[Size[i,PeriodCapture[i]]] 
    eReported[i,PeriodCapture[i]] <- bReported
    eReleased[i,PeriodCapture[i]] <- bReleased

    InLake[i,PeriodCapture[i]] <- 1
    Alive[i,PeriodCapture[i]] <- 1
    TBarTag[i,PeriodCapture[i]] ~ dbern(1-bTagLoss)

    Detected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eDetected[i,PeriodCapture[i]])
    DetectedCensus[i,PeriodCapture[i]] ~ dbern(InLake[i,PeriodCapture[i]] * Monitored[i,PeriodCapture[i]] * eDetectedCensus[i,PeriodCapture[i]] * Census[PeriodCapture[i]])
    MovedH[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMovedH[i,PeriodCapture[i]])
    MovedV[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Sensor[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMovedV[i,PeriodCapture[i]])
    
    Recaptured[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * eRecaptured[i,PeriodCapture[i]])
    Reported[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReported[i,PeriodCapture[i]] * TBarTag[i,PeriodCapture[i]])
    Released[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReleased[i,PeriodCapture[i]])

    for(a in (AnnualCapture[i]):nAnnual){
      eSpawning[i,a] <- (((LengthSpawned[i,a]^bSp) / (bLs^bSp + LengthSpawned[i,a]^bSp)) * bEs) 
    }
  for(j in (PeriodCapture[i]+1):nPeriod) {

    eMortalityHandling[i,j] <- ilogit(bMortalityHandling) * step(bMonthsHandling - eMonthsSinceCapture[i,j-1] - 1)
    eMonthsSinceCapture[i,j] <- (1-Recaptured[i,j]) * (eMonthsSinceCapture[i,j-1] + 1)
    eMortalitySpawning[i,j] <- ilogit(bMortalitySpawning) * Spawned[i,j]
    logit(eMortality[i,j]) <- bMortality + bMortalityAnnual[Annual[j]] + bMortalitySize[Size[i, j]] * (1 - Spawned[i,j]) + bMortalityMonth[Month[j]] * (1 - Spawned[i,j])
    Spawned[i,j] ~ dbern(Monitored[i,j] * SpawningPeriod[j] * eSpawning[i,Annual[j]])
    
    eDetected[i,j] <- Alive[i,j] * eDetectedAlive[i] + (1-Alive[i,j]) * eDetectedDead[i] 
    eDetectedCensus[i,j] <- bDetectedCensus

    eMovedH[i,j] <- bMovedH
    eMovedV[i,j] <- bMovedV
    logit(eRecaptured[i,j]) <-  bRecaptured + bRecapturedMonth[Month[j]] + bRecapturedSize[Size[i,j]]
    eReported[i,j] <- bReported
    eReleased[i,j] <- bReleased

    InLake[i,j] ~ dbern(InLake[i,j-1] * ((1 - (Recaptured[i,j-1])) + Recaptured[i,j-1] * Released[i,j-1]))
    Alive[i,j] ~ dbern(InLake[i,j] * Alive[i,j-1] * (1-eMortality[i,j-1]) * (1-eMortalityHandling[i,j-1]) * (1-eMortalitySpawning[i,j-1])) 
    TBarTag[i,j] ~ dbern(TBarTag[i,j-1] * (1-Recaptured[i,j-1]) * (1-bTagLoss))

    Detected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetected[i,j])
    DetectedCensus[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetectedCensus[i,j] * Census[j])
    MovedH[i,j] ~ dbern(Alive[i,j-1] * Monitored[i,j] * eMovedH[i,j])
    MovedV[i,j] ~ dbern(Alive[i,j] * Sensor[i,j] * Monitored[i,j] * eMovedV[i,j])
    
    Recaptured[i,j] ~ dbern(Alive[i,j] * eRecaptured[i,j])
    Reported[i,j] ~ dbern(Recaptured[i,j] * eReported[i,j] * TBarTag[i,j])
    Released[i,j] ~ dbern(Recaptured[i,j] * eReleased[i,j])
  }
  }

Block 1. Survival model description.

2024-04-23

Cowichan Lake Cuththroat Trout Exploitation Analysis 2022

Dalgarno, S. & Thorley, J.L. (2023) Cowichan Lake Cuththroat Trout Exploitation Analysis 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/86098464.

Background

The objectives of this analysis are to:

Estimate fishing and natural mortality.
Estimate annual variation in fishing and natural mortality.
Estimate seasonal variation in fishing and natural mortality.
Estimate effect of handling on natural mortality.
Estimate variation in fishing and natural mortality by fork length.

Data Preparation

The data were provided to Poisson Consulting Ltd. by Aswea Porter following extensive data screening.

Statistical Analysis

Model adequacy was assessed via posterior predictive checks (Kery and Schaub 2011). More specifically, the number of zeros and the first four central moments (mean, variance, skewness and kurtosis) for the deviance residuals were compared to the expected values by simulating new residuals. In this context the s-value indicates how surprising each observed metric is given the estimated posterior probability distribution for the residual variation.

The analyses were implemented using R version 4.2.2 (R Core Team 2021) and the mbr family of packages.

Model Descriptions

Survival

The natural mortality and recapture probability were estimated using a Bayesian individual state-space survival model (Thorley and Andrusak 2017) with monthly intervals. The survival model incorporated natural and handling mortality, acoustic detections and movements from a fixed receiver array and quarterly lake census, vertical movements from acoustic tags equipped with depth sensors, FLOY tag loss, recapture, reporting and release.

The expected fork length at monthly period $t$ was estimated from length at capture ($L_{i,fi}$) plus the length increment expected under a Von Bertalanffy Growth Curve (Walters and Martell 2004) \[L_{i,t} = L_{i,fi} + (L_{∞} − L_{i,fi} )(1 − e^{(−\frac{1}{12} \cdot k(t −fi))})\] where $fi$ is the period at capture and $\frac{1}{12}$ in the exponent adjusts for the fact that there are 12 periods per year.Linf was fixed to a biologically plausible value of 775 mm based on prior information (personal communication, R. Ptolemy); k of 0.23 was estimated from 88 aged Cutthroat Trout in Cowichan Lake (personal communication, Brendan Andersen).

In addition to assumptions 1, 4 to 10, and 13 in Thorley and Andrusak (2017), the survival model assumes that:

The effect of handling mortality is restricted to the first two monthly periods after initial capture or recapture.
Recapture probability varies by bi-monthly period and fork length.
Natural mortality varies by season, study year and size class.
The probability of detection depends on whether a fish is alive or has died near or far from a receiver and whether a census has occurred.
The monthly interval probability of FLOY tag loss is 0.001.
All recaptured fish that are released have their FLOY tags removed.
Reporting of a recaptured fish with a FLOY tag is likely greater than 90%.
Survival of fish tagged with a FLOY tag only is the same as those tagged with acoustic tags and a FLOY tag.
Acoustic tags are functional over their estimated battery lifespan and are not shed.

Seasons are defined as winter (December-February), spring (March-May), summer (June-August), and fall (September-November).

Study years span from October 01 to September 30 of the following year (i.e., study year 1 is 2019-10-01 to 2020-09-30). This was used instead of calendar year to ensure that each year group includes all months.

Size classes are defined as small ($\le$ 350 mm), medium (350 $-$ 500 mm) and large (> 500 mm).

Preliminary analyses indicated lack of support for variation in recapture probability by study year.

We assumed that the monthly interval probability of FLOY tag loss was 0.001 (0.012 annual probability) based on tag loss estimates from double-tagged Lake Trout in Horse Lake.

Model Templates

Survival

.model{
  bMortality ~ dnorm(-3, 3^-2)
  bMortalityHandling ~ dnorm(0, 2^-2)
  bMonthsHandling <- 2

  bDetectedAlive ~ dnorm(0, 2^-2)
  bDetectedAliveCensus ~ dnorm(0, 2^-2)
  bDieNear ~ dbeta(1, 1)
  bDetectedDeadNear ~ dnorm(1, 1^-2)
  bDetectedDeadNearCensus ~ dnorm(0, 2^-2)
  bDetectedDeadFar ~ dnorm(-1, 1^-2)
  bDetectedDeadFarCensus ~ dnorm(0, 2^-2)
  bMovedH ~ dbeta(1, 1)
  bMovedV ~ dbeta(1, 1)
  bRecaptured ~ dnorm(0, 2^-2)
  bReported ~ dbeta(1, 1)
  bReleased ~ dbeta(1, 1)
  # mode of 0.001
  bTagLoss <- 0.001
  bRecapturedLength ~ dnorm(0, 2^-2)

  sMortalitySizeSeason ~ dexp(1)
  for (i in 1:nSize){
    for(j in 1:nSeason){
      bMortalitySizeSeason[i,j] ~ dnorm(0, sMortalitySizeSeason^-2)
    }
  } 
  bMortalitySeason[1] <- 0
  for (i in 2:nSeason) {
    bMortalitySeason[i] ~ dnorm(0, 2^-2)
  }
  bMortalitySize[1] <- 0
  for (i in 2:nSize) {
    bMortalitySize[i] ~ dnorm(0, 2^-2)
  }
  bMortalityAnnual[1] <- 0
  for (i in 2:nAnnual) {
    bMortalityAnnual[i] ~ dnorm(0, 2^-2)
  }
  bRecapturedBiMonth[1] <- 0
  for (i in 2:nBiMonth) {
    bRecapturedBiMonth[i] ~ dnorm(0, 2^-2)
  }
  for (i in 1:nCapture){
    eMortalityHandling[i, PeriodCapture[i]] <- ilogit(bMortalityHandling)
    eMonthsSinceCapture[i,PeriodCapture[i]] <- 1-Recaptured[i,PeriodCapture[i]]
    
    logit(eDetectedAlive[i, PeriodCapture[i]]) <- bDetectedAlive 
    eDieNear[i] ~ dbern(bDieNear)
    logit(eDetectedDeadNear[i, PeriodCapture[i]]) <- bDetectedDeadNear 
    logit(eDetectedDeadFar[i, PeriodCapture[i]]) <- bDetectedDeadFar 
    eDetectedDead[i, PeriodCapture[i]] <- eDieNear[i] * eDetectedDeadNear[i, PeriodCapture[i]] + 
      (1 - eDieNear[i]) * eDetectedDeadFar[i, PeriodCapture[i]]
    
    logit(eMortality[i,PeriodCapture[i]]) <- bMortality + bMortalitySeason[Season[PeriodCapture[i]]] + bMortalityAnnual[Annual[PeriodCapture[i]]] + bMortalityHandling + bMortalitySize[Size[i, PeriodCapture[i]]] + bMortalitySizeSeason[Size[i, PeriodCapture[i]], Season[PeriodCapture[i]]]

    eMovedH[i,PeriodCapture[i]] <- bMovedH
    eMovedV[i,PeriodCapture[i]] <- bMovedV
    eDetected[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * eDetectedAlive[i,PeriodCapture[i]] +
      (1-Alive[i,PeriodCapture[i]]) * eDetectedDead[i,PeriodCapture[i]]
      
    logit(eRecaptured[i,PeriodCapture[i]]) <- bRecaptured + bRecapturedBiMonth[BiMonth[PeriodCapture[i]]] + bRecapturedLength * Length[i, PeriodCapture[i]] 
    eReported[i,PeriodCapture[i]] <- bReported
    eReleased[i,PeriodCapture[i]] <- bReleased

    InLake[i,PeriodCapture[i]] <- 1
    Alive[i,PeriodCapture[i]] ~ dbern((1-eMortality[i,PeriodCapture[i]]) * (1-eMortalityHandling[i,PeriodCapture[i]]))
    TBarTag[i,PeriodCapture[i]] ~ dbern(1-bTagLoss)

    Detected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eDetected[i,PeriodCapture[i]])
    MovedH[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMovedH[i,PeriodCapture[i]])
    MovedV[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Sensor[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMovedV[i,PeriodCapture[i]])
    
    Recaptured[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * eRecaptured[i,PeriodCapture[i]])
    Reported[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReported[i,PeriodCapture[i]] * TBarTag[i,PeriodCapture[i]])
    Released[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReleased[i,PeriodCapture[i]])

  for(j in (PeriodCapture[i]+1):nPeriod) {
    
    eMortalityHandling[i,j] <- ilogit(bMortalityHandling) * step(bMonthsHandling - eMonthsSinceCapture[i,j-1] - 1)
    eMonthsSinceCapture[i,j] <- (1-Recaptured[i,j]) * (eMonthsSinceCapture[i,j-1] + 1)
    
    logit(eDetectedAlive[i,j]) <- bDetectedAlive + bDetectedAliveCensus * Census[j]
    logit(eDetectedDeadNear[i,j]) <- bDetectedDeadNear + bDetectedDeadNearCensus * Census[j]
    logit(eDetectedDeadFar[i,j]) <- bDetectedDeadFar + bDetectedDeadFarCensus * Census[j]
    eDetectedDead[i,j] <- eDieNear[i] * eDetectedDeadNear[i,j] + (1 - eDieNear[i]) * eDetectedDeadFar[i,j]
    eDetected[i,j] <- Alive[i,j] * eDetectedAlive[i,j] + (1-Alive[i,j]) * eDetectedDead[i,j]

    logit(eMortality[i,j]) <- bMortality + bMortalityAnnual[Annual[j]] + bMortalitySize[Size[i,j]] + bMortalitySeason[Season[j]] + bMortalitySizeSeason[Size[i, j], Season[j]]

    eMovedH[i,j] <- bMovedH
    eMovedV[i,j] <- bMovedV
    logit(eRecaptured[i,j]) <-  bRecaptured + bRecapturedBiMonth[BiMonth[j]] + bRecapturedLength * Length[i,j]  
    eReported[i,j] <- bReported
    eReleased[i,j] <- bReleased

    InLake[i,j] ~ dbern(InLake[i,j-1] * ((1 - (Recaptured[i,j-1])) + Recaptured[i,j-1] * Released[i,j-1])) 
    Alive[i,j] ~ dbern(InLake[i,j] * Alive[i,j-1] * (1-eMortality[i,j]) * (1-eMortalityHandling[i,j])) 
    TBarTag[i,j] ~ dbern(TBarTag[i,j-1] * (1-Recaptured[i,j-1]) * (1-bTagLoss))

    Detected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetected[i,j])
    MovedH[i,j] ~ dbern(Alive[i,j-1] * Detected[i,j] * eMovedH[i,j])
    MovedV[i,j] ~ dbern(Alive[i,j] * Sensor[i,j] * Detected[i,j] * eMovedV[i,j])
    
    Recaptured[i,j] ~ dbern(Alive[i,j] * eRecaptured[i,j])
    Reported[i,j] ~ dbern(Recaptured[i,j] * eReported[i,j] * TBarTag[i,j])
    Released[i,j] ~ dbern(Recaptured[i,j] * eReleased[i,j])
  }
  }

Block 1. Survival model description.

2023-04-03

Elk River WCT Abundance Analysis 2022

Hill, N.E., Pearson, A.D., & Thorley, J.L. (2023) Elk River WCT Abundance Analysis 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/876409566.

Background

The Elk River is one of eight upper Kootenay River tributaries (Bull, Michel, Skookumchuck, St. Mary, upper Kootenay, White and Wigwam) that are designated as Class II due to the importance of their recreational fisheries for Westslope Cutthroat Trout (Oncorhynchus clarkii lewisi). The population abundance of Westslope Cutthroat Trout in the Elk River is uncertain. To improve population estimates a mark-recapture analysis using Passive Integrated Transponder (PIT) tags deployed by both boat electrofishing and guides was recommended (J. L. Thorley 2021).

After a successful pilot study in 2021, data collection via boat electrofishing and guides continued in 2022. The goal of the present analysis is to use the additional data to get an improved estimate of Westslope Cutthroat Trout population abundance in the Elk River.

Data Preparation

The data were provided by Nupqu Limited Partnership. The data were prepared for analysis using R version 4.2.2 (R Core Team 2022) and organized in a SQLite database. All river distances are based on the BC Freshwater Atlas layer.

Key assumptions of the data preparation included:

PIT Tag numbers were recorded correctly in instances where the data was not verifiable by the PIT Reader log.
PIT Reader log was taken to be correct in instances where the recorded PIT Tag code did not match the PIT Reader log.
2021 data were submitted in Mountain time; 2022 data were submitted in GMT-6.
Fork lengths were measured accurately.
GPS coordinates of outing start/end points were extended if they did not overlaps locations of fish captured on a given outing using the coordinates of a previous outing.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The estimates were produced using JAGS (Plummer 2003) and STAN (Carpenter et al. 2017). For additional information on Bayesian estimation the reader is referred to McElreath (2020).

Model adequacy was assessed via posterior predictive checks (Kery and Schaub 2011). More specifically, the proportion of zeros in the data and the first four central moments (mean, variance, skewness and kurtosis) in the deviance residuals were compared to the expected values by simulating new data based on the posterior distribution and assumed sampling distribution and calculating the deviance residuals.

Where computationally practical, the sensitivity of the posteriors to the choice of prior distributions was evaluated by increasing the standard deviations of all normal, half-normal and log-normal priors by an order of magnitude and then using $\hat{R}$ to evaluate whether the samples were drawn from the same posterior distribution (Joseph L. Thorley and Andrusak 2017).

The parameters are summarised in terms of the point estimate, lower and upper 95% compatibility limits (Rafi and Greenland 2020) and the surprisal s-value (Greenland 2019). The estimate is the median (50th percentile) of the MCMC samples while the 95% CLs are the 2.5th and 97.5th percentiles. The s-value indicates how surprising it would be to discover that the true value of the parameter is in the opposite direction to the estimate (Greenland 2019). An s-value of $>$ 4.32 bits, which is equivalent to a p-value $<$ 0.05 (Kery and Schaub 2011; Greenland and Poole 2013), indicates that the surprise would be equivalent to throwing at least 4.3 heads in a row.

Model selection was based the heuristic of directional certainty (Kery and Schaub 2011). Fixed effects were included if their s-value was $>$ 4.32 bits (Kery and Schaub 2011). Based on a similar argument, random effects were included if their standard deviation had lower 95% CLs $>$ 5% of the estimate.

The results are displayed graphically by plotting the modeled relationships between individual variables and the response 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 average values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). Primary explanatory variables are evaluated based on their estimated effect sizes with 95% CLs (Bradford, Korman, and Higgins 2005)

The analyses were implemented using R version 4.2.2 (R Core Team 2022) and the mbr family of packages.

Model Descriptions

Mark-Recapture

The data were analysed using a hierarchical Bayesian mark-recapture abundance model. Each zone was split into two kilometres sites except sites at the end of the zone were added to the previous site if the length was less than one kilometre. Captures from sites that had less then 10% of their length surveyed were not included in the analysis.

Key assumptions of the mark-recapture model include:

Lineal density varies randomly by year and zone within year.
There is no movement of marked fish between sites within a sampling season.
There is no mortality of fish within a sampling season.
The probability of capturing a marked or unmarked fish is the same.
All recaptured fish are correctly identified as being marked and there is no tag loss.
Capture efficiency varies by electrofishing effort.
Capture efficiency varies randomly by site visit.
The number of recaptures (of fish marked at that site in that year) and the total number of fish caught are each binomially distributed.

Model Templates

Mark-Recapture

.model{
  bDensity ~ dnorm(6, 2^-2)
  bEfficiency ~ dnorm(-5, 2^-2)
  bEfficiencySecondsPerMetre ~ dnorm(0, 1^-2)
  bTheta ~ dexp(10)
  bDensityAnnual[1] <- 0
  for (i in 2:nannual) {
    bDensityAnnual[i] ~ dnorm(0, 2^-2)
  }
  sDensityAnnualZoneID ~ dnorm(0, 1^-2) T(0,)
  for (i in 1:nannual) {
    for (j in 1:nzone_id) {
      bDensityAnnualZoneID[i, j] ~ dnorm(0, sDensityAnnualZoneID^-2)
    }
  }
  bEfficiencyAnnual[1] <- 0
  for (i in 2:nannual) {
    bEfficiencyAnnual[i] ~ dnorm(0, 2^-2)
  }
  for (i in 1:nObs) {
    log(eDensity[i]) <- bDensity + bDensityAnnual[annual[i]] + bDensityAnnualZoneID[annual[i], zone_id[i]]
    logit(eProb[i]) <- bEfficiency + bEfficiencySecondsPerMetre * seconds_per_metre[i] + bEfficiencyAnnual[annual[i]]
    eEfficiency[i] ~ dbeta((2 * eProb[i]) / bTheta, (2 * (1 - eProb[i])) / bTheta)  
    recaptures[i] ~ dbin(eEfficiency[i], marked[i])
    bAbundance[i] ~ dpois(eDensity[i] * site_length[i] * survey_proportion[i])
    count[i] ~ dbin(eEfficiency[i], bAbundance[i])
  }

Block 1. Model description.

2023-03-30