The suggested citation for this analytic appendix is:
Amies-Galonski, E., Andrusak, G.F., Stephens, C.R.A., and Thorley, J.L. (2024) Kaslo Bull Trout Productivity 2023. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1402985678.
The Kaslo River and Keen Creek, which is a tributary of the Kaslo River, are important Bull Trout spawning and rearing tributaries on Kootenay Lake. From 2012 to 2023, field crews have night-snorkeled these systems in the fall and recorded all Bull Trout less than 350 mm in length. Keen Creek has not been snorkelled since 2019 due to visibility issues. Snorkel and electrofishing marking crews have also captured and tagged juvenile Bull Trout for the snorkel crews to resight. Redd counts have been conducted in both systems since 2006, with the exception of 2020, when Keen Creek was not surveyed. The primary goal of the current analyses is to answer the following questions:
What is the observer efficiency when night-snorkeling for juvenile Bull Trout in the Kaslo River and Keen Creek?
What are the numbers of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the stock (number of redds and eggs) and the resultant numbers of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 4.4.1 (R Core Team 2024).
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 et al. 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 summarized in terms of the point estimate, lower and upper 95% compatibility limits (Rafi and Greenland 2020) and the surprisal s-value (Greenland 2019). The coefficient of variation is also shown for juvenile abundance to understand sampling efforts. 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.
Variable selection was based on 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 median estimate.
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.
Where computationally practical, the sensitivity of the parameters 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 (Thorley and Andrusak 2017).
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). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford et al. 2005).
The analyses were implemented using R version 4.4.1
(R Core Team 2024) and the
mbr family of packages.
The annual bias (inaccuracy) and error (imprecision) in observer’s fish length estimates when spotlighting (standing) and snorkeling were quantified from the divergence of their length distribution from the length distribution for JLT and SH (the two most experience snorkelers) in that year. More specifically, the length correction that minimized the Jensen-Shannon divergence (Lin 1991) provided a measure of the inaccuracy while the minimum divergence (the Jensen-Shannon divergence was calculated with log to base 2 which means it lies between 0 and 1) provided a measure of the imprecision.
All resighted fish with a tag were allocated to the closest unallocated marked fish (with the same colour tag) by fork length and distance. The marked fish were analysed using a Bayesian logistic regression model. The key assumption of the logistic regression model is that:
The preliminary analysis for observer efficiency indicated that system, observer, gradient, sinuosity, and river kilometre were not informative predictors.
Both systems were broken into 100 m sites by bank. The lineal density at each site was estimated using an over-dispersed Poisson Generalized Linear Mixed Model.
Key assumptions of the Bayesian GLMM include:
The preliminary analysis for density indicated that site sinuosity and gradient were not informative predictors of lineal density.
The condition of adults was estimated using an allometric weight-length relationship. Key assumptions of the condition model include:
Female spawner fecundity was estimated using an allometric egg-weight relationship. The key assumptions of the fecundity model include:
The average length of spawners in each system and year was estimated using a linear regression. Key assumptions of the spawner length model include:
The total egg deposition in each year was estimated by
The stock-recruitment relationship between the total number eggs and the abundance of age-1 individuals was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
The \(E_{K/2}\) Limit Reference Point (Mace 1994) (\(E_{0.5 R_{max}}\)) was calculated, corresponding to the stock (number of eggs per 100 meters) that produce 50% of the maximum recruitment (\(K\)).
.model{
bIntercept ~ dnorm(0, 2^-2)
bLength ~ dnorm(0, 2^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Block 1. Final model.
.model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, 1)
b0 ~ dnorm(0, 5^-2)
bRkm ~ dnorm(0, 2^-2)
bSystem[1] <- 0
for(i in 2:nSystem) {
bSystem[i] ~ dnorm(0, 2^2)
}
sSystemYear ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, sSystemYear^-2)
}
}
sSite ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
sDispersion ~ dnorm(0, 2^-2) T(0,)
for(i in 1:length(Count)) {
log(eDensity[i]) <- b0 + bSystem[System[i]] + bRkm * Rkm[i] + bSite[Site[i]] + bSystemYear[System[i],Year[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Block 2. The final model.
.model{
b0 ~ dnorm(6, 2^-2)
bLength ~ dnorm(3, 1^-2)
sWeight ~ dnorm(0, 1^-2) T(0,)
sYear ~ dnorm(0, 1^-2) T(0,)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nObs) {
log(eWeight[i]) <- b0 + (log(Length[i]) - log(500)) * bLength + bYear[Year[i]]
Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
}Block 3. Model description.
.model{
b0 ~ dnorm(0, 2^-2)
bWeight ~ dnorm(1, 2^-2)
sFecundity ~ dnorm(0, 2^-2) T(0,)
for (i in 1:nObs) {
log(eFecundity[i]) <- b0 + (log(Weight[i]) - log(2300)) * bWeight
Fecundity[i] ~ dlnorm(log(eFecundity[i]), sFecundity^-2)
}Block 4. Model description.
.model{
bLength ~ dnorm(650, 100^-2)
sLength ~ dnorm(0, 100^-2) T(0,)
sYear ~ dnorm(0, 100^-2) T(0,)
sSystem ~ dnorm(0, 100^-2) T(0,)
sYearSystem ~ dnorm(0, 100^-2) T(0,)
bSex ~ dbeta(1, 1)
bFemale ~ dnorm(0, 100^-2)
bBias ~ dnorm(0, 100^-2)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nSystem) {
bSystem[i] ~ dnorm(0, sSystem^-2)
}
for (i in 1:nYear) {
for (j in 1:nSystem) {
bYearSystem[i, j] ~ dnorm(0, sYearSystem^-2)
}
}
bCatchType[1] <- 0
for(i in 2:nCatchType) {
bCatchType[i] ~ dnorm(0, 100^-2)
}
for (i in 1:nObs) {
Female[i] ~ dbern(bSex)
eLength[i] <- bLength + bSystem[System[i]] + bFemale * Female[i] + bCatchType[CatchType[i]] + bYear[Year[i]] + bYearSystem[Year[i], System[i]] + bBias * Bias[i]
Length[i] ~ dnorm(eLength[i], sLength^-2)
}Block 5. Model description.
.model {
bAlpha ~ dbeta(4, 49)
bK ~ dnorm(3, 1^-2)
bKPopulation ~ dnorm(0, 1^-2)
sRecruits ~ dexp(1)
for(i in 1:nObs){
log(eK[i]) <- bK + (Kaslo[i] * bKPopulation - (1-Kaslo[i]) * bKPopulation)
eBeta[i] <- bAlpha/eK[i]
eRecruits[i] <- bAlpha * Stock[i] / (1 + eBeta[i] * Stock[i])
Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
}Block 6. Final model.
The suggested citation for this analytic appendix is:
Thorley, J.L. (2024) Lower Columbia River Fish Population Indexing 2023. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/25988200.
In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.
The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:
What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.
The suggested citation for this analytic appendix is:
Thorley, J.L., Dalgarno, S., Lyons, S. (2023) Lower Columbia River Fish Population Indexing 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1392245070.
In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.
The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:
What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.
The suggested citation for this analytic appendix is:
Dalgarno, S. & Thorley, J.L. (2023) Quesnel Exploitation Analysis 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1230088647.
Quesnel Lake supports a recreational fishery for Rainbow Trout. To provide information on the natural and fishing mortality, Rainbow Trout were caught by angling and tagged with acoustic transmitters and/or a $100 and $10 reward tag.
The objectives of this analysis are to:
The outing, receiver deployment, detection, fish capture and recapture information were provided by the Ministry of Forests, Lands and Natural Resource Operations and added to a SQLite database.
The data were prepared for analysis using R version 4.3.0 (R Core Team 2022).
Receivers were assumed to have a detection range of 500 m. Detections were aggregated daily, where for each transmitter the receiver with the most number of detections was chosen. In the case of a tie, the receiver with the greatest coverage area was chosen. Only individuals with a fork length (FL) \(\geq\) 450 mm, an acoustic tag life \(\geq\) 365 days (if acoustically tagged) and a $100 and $10 reward tags were included in the survival 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)
Where possible, 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 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.0
(R Core Team 2022) and the
mbr family of packages.
The expected weight of fish of a given length were estimated from the data using a mass-length model (He et al. 2008).
More specifically the model was based on the allometric relationship
\[ W = \alpha L^{\beta}\]
where \(W\) is the weight (mass), \(\alpha\) is the coefficent, \(\beta\) is the exponent and \(L\) is the length.
To improve chain mixing the relation was log-transformed, i.e.,
\[ \log(W) = \log(\alpha) + \beta \cdot \log(L)\] Length and weights of fish were taken from two sources: captures from broodstock collection in the Horsefly River (2015-2022) and in-lake captures from tagging outings (2013-2022).
Key assumptions of the condition model include:
Annual growth of fish were estimated from inter-annual recaptures (in-lake and broodstock) and aged broodstock captures using the Fabens method (Fabens 1965) for estimating the von Bertalanffy growth curve (von Bertalanffy 1938). This curve is based on the premise that:
\[ \frac{\text{d}L}{\text{d}t} = k (L_{\infty} - L)\]
where \(L\) is the length of the individual, \(k\) is the growth coefficient and \(L_{\infty}\) is the maximum length.
Integrating the above equation gives:
\[ L_t = L_{\infty} (1 - e^{-k(t - t_0)})\]
where \(L_t\) is the length at time \(t\) and \(t_0\) is the time at which the individual would have had zero length.
The Fabens form allows
\[ L_r = L_c + (L_{\infty} - L_c) (1 - e^{-kT})\]
where \(L_r\) is the length at recapture, \(L_c\) is the length at capture and \(T\) is the time between capture and recapture.
For aged fish, length at capture is the length at Age-0 \(L_0\), estimated
with the modified form:
\[L_r = L_{0} - (L_{\infty} - L_0) (1 - e^{-kA})\] where \(A\) is the age.
Key assumptions of the growth model include:
A growth year spans April 01 to March 31 the following year.
Preliminary analyses indicated lack of support for variation in \(k\) by Horsefly River sockeye run size with a one or two year lag.
The natural mortality and recapture probability were estimated using a Bayesian individual state-space survival model (Thorley and Andrusak 2017) with three month intervals. The survival model incorporated natural and handling mortality, acoustic detections, inter-section movement, T-bar tag loss, spawning, recapture, reporting, fork length at spawning and sockeye run size in the previous year.
T-bar tag loss was estimated from the number of tags reported from recaught double-tagged individuals (Fabrizio et al. 1999). The model assumes that tag loss is not a time-dependent process and occurs following initial release.
Fork length at spawning was calculated from the measured length at capture (\(L_{i,fi}\)) plus the length increment expected under a Von Bertalanffy Growth Curve (Walters & Martell, 2004)
\[L_{i,t} = L_{i,fi} + (L_{∞} − L_{i,fi} )(1 − e^{(−0.25 \cdot k(t −fi))})\]
where the 0.25 in the exponent adjusts for the fact that there are four time periods per year. Values for \(L_\infty\) and annual \(k\) and were estimated in the growth model.
The probability of spawning at estimated length \(L\) is determined by
\[S = \frac{L^{S_p}}{L_s^{S_p} + L^{S_p}} \cdot es\]
where \(L_s\) is the length at which 50% of fish are mature, \(es\) is the probability of a mature fish spawning and \(S_p\) is the maturity (as a function of length) power.
Key assumptions of the spawning probability sub-model include:
A fish was considered to have spawned in a given year if one of the following conditions were met: detected within a spawning river during the spawning period with a detection hiatus (no detections in lake) of at least 7 days; or detected within a spawning gate during the spawning period with a detection hiatus of at least 14 days. Cutoff values for minimum hiatus periods were determined by inspecting the distribution of in-lake hiatus periods for fish detected in spawning rivers and the distribution of all hiatus periods for suspected spawners within the spawning period.
In addition to assumptions 1 to 2 and 4 to 10, in Thorley and Andrusak (2017), the model also assumes that:
Seasons are defined as winter (January to March), spring (April to June), summer (July to September), and fall (October to December).
Size classes are defined as small (\(\le\) 600 mm), medium (601 \(-\) 700 mm) and large (\(>\) 700 mm).
Sockeye abundance is the Horsefly River run size with a two year lag (i.e., spawning in the spring of 2016 is associated with autumn 2014 sockeye abundance).
Preliminary analyses indicated lack of support for variation in mortality by year, sockeye year (October 01 to September 30 the following year) and size class within year.
The in-lake and spawner abundance for individuals \(\ge\) 500 mm were estimated using hierarchical Bayesian mark-recapture sub-models within the survival model, where the number of alive marked individuals in the lake at each time period was estimated from the ‘Alive’ latent state in the previous period. The number of spawning marked individuals was estimated from the spawning probability and the number of alive marked individuals in the previous period.
The ‘in-lake’ mark-recapture model estimated abundance from in-lake captures and recaptures from researcher tagging outings.
The ‘spawner’ mark-recapture model estimated spawner abundance from spawner captures and recaptures from broodstock collection in the Horsefly River.
An additional in-lake abundance \(A_L\) estimate was derived from estimated spawner abundance \(A_S\) and spawning probability \(P_S\) given the following: \[A_L = \frac{A_S}{P_S} \]
Key assumptions of the mark-recapture models include:
Additional assumptions for the lake model include:
The optimal yield (defined as the harvested number of fish as percent of the number of recruits at carrying capacity) was calculated using a yield-per-recruit approach (Bison, O’Brien, and Martell 2003).
Scenarios explored include:
Llo) of 40, 50, 60, 70, 80 and 90, with other
YPR parameter values held constant.Rk (lifetime spawners per spawner at low density) values of 5,
7.5, and 10 and Ls (length at 50% mature) values of 55, 60, and
65, with Llo held at 60.Rk values of 5, 7.5, and 10 and Ls values of 55, 60, and 65,
with Llo held at 50.Key assumptions include:
Llo are retained.Llo are released..model{
bWeight ~ dnorm(3, 2^-2)
bWeightLength ~ dnorm(3, 1^-2)
bWeightDayte ~ dnorm(0, 1^-2);
bWeightLengthDayte ~ dnorm(0, 1^-2)
sWeightAnnual ~ dexp(1)
for(i in 1:nAnnual) {
bWeightAnnual[i] ~ dnorm(0, sWeightAnnual^-2)
}
sWeightSourceAnnual ~ dexp(1)
for(i in 1:nSource){
for(j in 1:nAnnual){
bWeightSourceAnnual[i,j] ~ dnorm(0, sWeightSourceAnnual^-2)
}
}
sWeight ~ dexp(1)
for(i in 1:nObs) {
log(eWeight[i]) <- bWeight + bWeightDayte * Dayte[i] + bWeightAnnual[Annual[i]] + bWeightSourceAnnual[Source[i], Annual[i]] + (bWeightLength + bWeightLengthDayte * Dayte[i]) * (log(Length[i]) - log(100))
Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
}Block 1. Model description.
.model {
bK ~ dnorm(0, 2^-2)
bL0 ~ dnorm(0, 100^-2)
sKYear ~ dexp(1)
theta <- 0.5
bLinf ~ dnorm(900, 100^-2) T(0,)
sGrowth ~ dexp(1)
sLength ~ dexp(1)
for (i in 1:nYear) {
bKYear[i] ~ dnorm(0, sKYear^-2)
log(eK[i]) <- bK + bKYear[i]
}
for(i in 1:nAged){
eLength[i] <- bLinf - (bLinf - bL0) * exp(-exp(bK) * Age[i])
Length[i] ~ dnorm(eLength[i], sLength^-2)
}
for(i in 1:nFish){
eL0[i] <- max(bL0, LengthAtRelease[i])
for(j in 1:nYear){
ekProp[i,j] <- eK[j] * pYear[i,j]
}
eSumEk[i] <- exp(-sum(ekProp[i,]))
eGrowth[i] <- max(0, (bLinf - eL0[i]) * (1 -eSumEk[i]))
Growth[i] ~ dt(eGrowth[i], sGrowth^-2, 1/theta)
}Block 2.
model {
bTagLoss ~ dbeta(1, 1)
for (i in 1:nObs) {
TagsRecap[i] ~ dbin(1 - bTagLoss, 2) T(1,)
}Block 3. Model description.
.model{
bMortality ~ dnorm(-3, 3^-2)
bMortalitySpawning ~ dnorm(0, 2^-2)
bLs ~ dnorm(6, 2^-2)
bSp ~ dexp(1/20)
bEs ~ dnorm(0, 2^-2)
bDetectedAlive ~ dbeta(1, 1)
bDieNear ~ dbeta(1, 1)
bDetectedDeadNear ~ dbeta(10, 1)
bDetectedDeadFar ~ dbeta(1, 10)
bMoved ~ dbeta(1, 1)
bRecaptured ~ dnorm(0, 2^-2)
bReported ~ dbeta(10, 1)
bReleased ~ dbeta(1, 1)
bReleasedResearcher ~ dbeta(1, 1)
bRecapturedResearcher ~ dbeta(5, 1)
bTagLoss <- 0.0545
bTagLossT <- 0.0001
bEfficiency ~ dnorm(-4, 2^-2)
bAbundance ~ dnorm(7, 2^-2)
bOutingHours ~ dnorm(0, 2^-2)
bEfficiencySpawn ~ dnorm(-4, 2^-2)
bAbundanceSpawn ~ dnorm(1000, 500^-2)
bMortalitySeason[1] <- 0
for (i in 2:nSeason) {
bMortalitySeason[i] ~ dnorm(0, 2^-2)
}
for (i in 1:nSizeClass) {
bMortalitySockeye[i] ~ dnorm(0, 2^-2)
}
sRecapturedAnnual ~ dexp(1)
for (i in 1:nAnnual) {
bRecapturedAnnual[i] ~ dnorm(0, sRecapturedAnnual^-2)
}
bRecapturedSeason[1] <- 0
for (i in 2:nSeason) {
bRecapturedSeason[i] ~ dnorm(0, 2^-2)
}
sLsAnnual ~ dexp(1)
for (i in 1:nAnnual) {
bLsAnnual[i] ~ dnorm(0, sLsAnnual^-2)
}
sEsAnnual ~ dexp(1)
for (i in 1:nAnnual) {
bEsAnnual[i] ~ dnorm(0, sEsAnnual^-2)
}
for (i in 1:nCapture){
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 + bMortalitySeason[Season[PeriodCapture[i]]] + bMortalitySockeye[SizeClass[i,PeriodCapture[i]]] * SockeyePrevious[PeriodCapture[i]]
logit(eRecaptured[i,PeriodCapture[i]]) <- bRecaptured + bRecapturedSeason[Season[PeriodCapture[i]]] + bRecapturedAnnual[Annual[PeriodCapture[i]]]
eMoved[i,PeriodCapture[i]] <- bMoved
eReported[i,PeriodCapture[i]] <- bReported
eReleased[i,PeriodCapture[i]] <- bReleased
eReleasedResearcher[i,PeriodCapture[i]] <- bReleasedResearcher
eRecapturedResearcher[i,PeriodCapture[i]] <- bRecapturedResearcher
InLake[i,PeriodCapture[i]] <- 1
Alive[i,PeriodCapture[i]] ~ dbern(1-eMortality[i,PeriodCapture[i]])
Alive50[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * (Length[i,PeriodCapture[i]] >= 500)
TBarTag100[i,PeriodCapture[i]] ~ dbern(1-bTagLoss)
TBarTag10[i,PeriodCapture[i]] ~ dbern(1-bTagLoss)
eDetected[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * eDetectedAlive[i] + (1-Alive[i,PeriodCapture[i]]) * eDetectedDead[i]
Detected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eDetected[i,PeriodCapture[i]])
Moved[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMoved[i,PeriodCapture[i]])
RecapturedResearcher[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * eRecapturedResearcher[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]] * step(TBarTag100[i,PeriodCapture[i]] + TBarTag10[i,PeriodCapture[i]] - 1))
Released[i,PeriodCapture[i]] ~ dbern(Recaptured[i,PeriodCapture[i]] * eReleased[i,PeriodCapture[i]])
ReleasedResearcher[i,PeriodCapture[i]] ~ dbern(RecapturedResearcher[i,PeriodCapture[i]] * eReleasedResearcher[i,PeriodCapture[i]])
for(j in (PeriodCapture[i]+1):nPeriod) {
log(eLs[i,j]) <- bLs + bLsAnnual[Annual[j]]
logit(eEs[i,j]) <- bEs + bEsAnnual[Annual[j]]
eSpawning[i,j] <- (((Length[i,j]^bSp) / (eLs[i,j]^bSp + Length[i,j]^bSp)) * eEs[i,j])
Spawned[i,j] ~ dbern(Alive[i,j-1] * Monitored[i,j] * SpawningPeriod[j] * eSpawning[i,j])
logit(eMortality[i,j]) <- bMortality + bMortalitySpawning * Spawned[i,j-1] + bMortalitySeason[Season[j]] + bMortalitySockeye[SizeClass[i,j]] * SockeyePrevious[j]
logit(eRecaptured[i,j]) <- bRecaptured + bRecapturedSeason[Season[j]] + bRecapturedAnnual[Annual[j]]
eMoved[i,j] <- bMoved
eRecapturedResearcher[i,j] <- bRecapturedResearcher
eReported[i,j] <- bReported
eReleased[i,j] <- bReleased
eReleasedResearcher[i,j] <- bReleasedResearcher
InLake[i,j] ~ dbern(InLake[i,j-1] * ((1 - (Recaptured[i,j-1])) + Recaptured[i,j-1] * Released[i,j-1]) * ((1 - (RecapturedResearcher[i,j-1])) + RecapturedResearcher[i,j-1] * ReleasedResearcher[i,j-1]))
Alive[i,j] ~ dbern(InLake[i,j] * Alive[i,j-1] * (1-eMortality[i,j]))
Alive50[i,j] <- Alive[i,j] * (Length[i,j] >= 500)
TBarTag100[i,j] ~ dbern(TBarTag100[i,j-1] * (1-Recaptured[i,j-1]) * (1-bTagLossT))
TBarTag10[i,j] ~ dbern(TBarTag10[i,j-1] * (1-Recaptured[i,j-1]) * (1-bTagLossT))
eDetected[i,j] <- Alive[i,j] * eDetectedAlive[i] + (1-Alive[i,j]) * eDetectedDead[i]
Detected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetected[i,j])
Moved[i,j] ~ dbern(Alive[i,j] * Detected[i,j] * eMoved[i,j])
RecapturedResearcher[i,j] ~ dbern(Alive[i,j] * eRecapturedResearcher[i,j] * step(TBarTag100[i,j] + TBarTag10[i,j] - 1))
Recaptured[i,j] ~ dbern(Alive[i,j] * eRecaptured[i,j])
Reported[i,j] ~ dbern(Recaptured[i,j] * eReported[i,j] * step(TBarTag100[i,j] + TBarTag10[i,j] - 1))
Released[i,j] ~ dbern(Recaptured[i,j] * eReleased[i,j])
ReleasedResearcher[i,j] ~ dbern(RecapturedResearcher[i,j] * eReleasedResearcher[i,j])
}
}
Marked[1] <- sum(Alive50[1:PeriodCaptures[1],1])
MarkedSpawners[1] <- 1
for(p in 2:nPeriod){
Marked[p] <- sum(Alive50[1:PeriodCaptures[p],p])
logit(eSpawningPeriod[p]) <- bEs + bEsAnnual[Annual[p]]
MarkedSpawners[p] ~ dpois(Marked[p-1] * eSpawningPeriod[p])
}
sAbundanceAnnual ~ dexp(1)
for (i in 1:nAnnualOuting) {
bAbundanceAnnual[i] ~ dnorm(0, sAbundanceAnnual^-2)
}
eAbundanceSpawn ~ dpois(bAbundanceSpawn)
for(l in 1:nPeriodCount){
logit(eEfficiencySpawn[l]) <- bEfficiencySpawn
RecaptureSpawners[l] ~ dbin(eEfficiencySpawn[l], MarkedSpawners[PeriodCount[l]])
CaptureSpawners[l] ~ dbin(eEfficiencySpawn[l], eAbundanceSpawn)
}
for(k in 1:nPeriodOuting){
log(eAbundance[k]) <- bAbundance + bAbundanceAnnual[AnnualOuting[k]]
eAbundanceInt[k] ~ dpois(eAbundance[k])
logit(eEfficiency[k]) <- bEfficiency + bOutingHours * OutingHours[PeriodOuting[k]]
Recaptures[k] ~ dbin(eEfficiency[k], Marked[PeriodOuting[k]-1])
Captures[k] ~ dbin(eEfficiency[k], eAbundanceInt[k])
}Block 4. Survival model description.
The suggested citation for this analytic appendix is:
Amies-Galonski, E., Andrusak G.F. and Thorley, J.L. (2023) Kaslo Bull Trout Productivity 2022. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/798364147.
The Kaslo River and Keen Creek, which is a tributary of the Kaslo River, are important Bull Trout spawning and rearing tributaries on Kootenay Lake. From 2012 to 2022, field crews have night-snorkeled these systems in the fall and recorded all Bull Trout less than 350 mm in length. Keen Creek was not snorkelled in 2020, 2021, or 2022 due to visibility issues. Snorkel and electrofishing marking crews have also captured and tagged juvenile Bull Trout for the snorkel crews to resight. Redd counts have been conducted in both systems since 2006, with the exception of 2020, when Keen Creek was not surveyed. The primary goal of the current analyses is to answer the following questions:
What is the observer efficiency when night-snorkeling for juvenile Bull Trout in the Kaslo River and Keen Creek?
What are the numbers of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the stock (number of redds and eggs) and the resultant numbers of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 4.3.0 (R Core Team 2022).
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 et al. 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% 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.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.
Variable selection was based on 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 median estimate.
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.
Where computationally practical, the sensitivity of the parameters 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 (Thorley and Andrusak 2017).
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). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford et al. 2005).
The analyses were implemented using R version 4.3.0
(R Core Team 2022) and the
mbr family of packages.
The annual bias (inaccuracy) and error (imprecision) in observer’s fish length estimates when spotlighting (standing) and snorkeling were quantified from the divergence of their length distribution from the length distribution for JLT and SH (the two most experience snorkelers) in that year. More specifically, the length correction that minimised the Jensen-Shannon divergence (Lin 1991) provided a measure of the inaccuracy while the minimum divergence (the Jensen-Shannon divergence was calculated with log to base 2 which means it lies between 0 and 1) provided a measure of the imprecision.
All resighted fish with a tag were allocated to the closest unallocated marked fish (with the same colour tag) by fork length and distance. The marked fish were analysed using a Bayesian logistic regression model. The key assumption of the logistic regression model is that:
The preliminary analysis for observer efficiency indicated that system, observer, gradient, sinuosity, and river kilometre were not informative predictors.
Both systems were broken into 100 m sites by bank. The lineal density at each site was estimated using an over-dispersed Poisson Generalized Linear Mixed Model.
Key assumptions of the Bayesian GLMM include:
The preliminary analysis for density indicated that site sinuosity and gradient were not informative predictors of lineal density.
The condition of adults was estimated using an allometric weight-length relationship. Key assumptions of the condition model include:
Female spawner fecundity was estimated using an allometric egg-weight relationship. The key assumptions of the fecundity model include:
The average length of spawners in each system and year was estimated using a linear regression. Key assumptions of the spawner length model include:
The total egg deposition in each year was estimated by
The stock-recruitment relationship between the total number eggs and the abundance of age-1 individuals was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
The \(E_{K/2}\) Limit Reference Point (Mace 1994) (\(E_{0.5 R_{max}}\)) was calculated, corresponding to the stock (number of eggs per 100 meters) that produce 50% of the maximum recruitment (\(K\)).
model{
bIntercept ~ dnorm(0, 2^-2)
bLength ~ dnorm(0, 2^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Block 1. Final model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, 1)
b0 ~ dnorm(0, 5^-2)
bRkm ~ dnorm(0, 2^-2)
bSystem[1] <- 0
for(i in 2:nSystem) {
bSystem[i] ~ dnorm(0, 2^2)
}
sSystemYear ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, sSystemYear^-2)
}
}
sSite ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
sDispersion ~ dnorm(0, 2^-2) T(0,)
for(i in 1:length(Count)) {
log(eDensity[i]) <- b0 + bSystem[System[i]] + bRkm * Rkm[i] + bSite[Site[i]] + bSystemYear[System[i],Year[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Block 2. The final model.
model{
b0 ~ dnorm(6, 2^-2)
bLength ~ dnorm(3, 1^-2)
sWeight ~ dnorm(0, 1^-2) T(0,)
sYear ~ dnorm(0, 1^-2) T(0,)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nObs) {
log(eWeight[i]) <- b0 + (log(Length[i]) - log(500)) * bLength + bYear[Year[i]]
Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
}Block 3. Model description.
model{
b0 ~ dnorm(0, 2^-2)
bWeight ~ dnorm(1, 2^-2)
sFecundity ~ dnorm(0, 2^-2) T(0,)
for (i in 1:nObs) {
log(eFecundity[i]) <- b0 + (log(Weight[i]) - log(2300)) * bWeight
Fecundity[i] ~ dlnorm(log(eFecundity[i]), sFecundity^-2)
}Block 4. Model description.
model{
bLength ~ dnorm(650, 100^-2)
sLength ~ dnorm(0, 100^-2) T(0,)
sYear ~ dnorm(0, 100^-2) T(0,)
sSystem ~ dnorm(0, 100^-2) T(0,)
sYearSystem ~ dnorm(0, 100^-2) T(0,)
bSex ~ dbeta(1, 1)
bFemale ~ dnorm(0, 100^-2)
bBias ~ dnorm(0, 100^-2)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nSystem) {
bSystem[i] ~ dnorm(0, sSystem^-2)
}
for (i in 1:nYear) {
for (j in 1:nSystem) {
bYearSystem[i, j] ~ dnorm(0, sYearSystem^-2)
}
}
bCatchType[1] <- 0
for(i in 2:nCatchType) {
bCatchType[i] ~ dnorm(0, 100^-2)
}
for (i in 1:nObs) {
Female[i] ~ dbern(bSex)
eLength[i] <- bLength + bSystem[System[i]] + bFemale * Female[i] + bCatchType[CatchType[i]] + bYear[Year[i]] + bYearSystem[Year[i], System[i]] + bBias * Bias[i]
Length[i] ~ dnorm(eLength[i], sLength^-2)
}Block 5. Model description.
model {
bAlpha ~ dbeta(4, 49)
bK ~ dnorm(3, 1^-2)
bKPopulation ~ dnorm(0, 1^-2)
sRecruits ~ dexp(1)
for(i in 1:nObs){
log(eK[i]) <- bK + (Kaslo[i] * bKPopulation - (1-Kaslo[i]) * bKPopulation)
eBeta[i] <- bAlpha/eK[i]
eRecruits[i] <- bAlpha * Stock[i] / (1 + eBeta[i] * Stock[i])
Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
}Block 6. Final model.
The suggested citation for this analytic appendix is:
Dalgarno, S. & Thorley, J.L. (2022) Quesnel Exploitation Analysis 2021. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/978976227.
Quesnel Lake supports a recreational fishery for Rainbow Trout. To provide information on the natural and fishing mortality, Rainbow Trout were caught by angling and tagged with acoustic transmitters and/or a $100 and $10 reward tag.
The outing, receiver deployment, detection, fish capture and recapture information were provided by the Ministry of Forests, Lands and Natural Resource Operations and added to a SQLite database.
The data were prepared for analysis using R version 4.2.1 (R Core Team 2021).
Receivers were assumed to have a detection range of 500 m. Detections were aggregated daily, where for each transmitter the receiver with the most number of detections was chosen. In the case of a tie, the receiver with the greatest coverage area was chosen. Only individuals with a fork length (FL) \(\geq\) 450 mm, an acoustic tag life \(\geq\) 365 days (if acoustically tagged) and a $100 and $10 reward tags were included in the survival 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 (2016).
Unless stated otherwise, the Bayesian analyses used weakly informative normal and half-normal or uninformative beta 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 can be considered a test of directionality. More specifically it indicates how surprising (in bits) it would be to discover that the true value of the parameter is in the opposite direction to the estimate. An s-value of 4.3 bits, which is equivalent to a p-value (Kery and Schaub 2011; Greenland and Poole 2013) of 0.05, indicates that the surprise would be equivalent to throwing 4.3 heads in a row. The condition that non-essential explanatory variables have s-values \(\geq\) 4.3 bits provides a useful model selection heuristic (Kery and Schaub 2011).
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). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent or n-fold change in the response variable) with 95% credible intervals (CIs, Bradford, Korman, and Higgins 2005).
The analyses were implemented using R version 4.2.1
(R Core Team 2020) and the
mbr family of packages.
The expected weight of fish of a given length were estimated from the data using a mass-length model (He et al. 2008).
More specifically the model was based on the allometric relationship
\[ W = \alpha L^{\beta}\]
where \(W\) is the weight (mass), \(\alpha\) is the coefficent, \(\beta\) is the exponent and \(L\) is the length.
To improve chain mixing the relation was log-transformed, i.e.,
\[ \log(W) = \log(\alpha) + \beta \cdot \log(L)\]
Key assumptions of the condition model include:
The expected length of fish at a given age were estimated from broodstock data collected in the Horsefly River from 2015-2021 using a Von Bertalanffy growth curve model (von Bertalanffy 1938).
\[ L = L_\infty \cdot (1 - \exp(-k \cdot (A - t_0) ))\]
where \(A\) is the scale age in years and \(t_0\) is the extrapolated age at zero length.
Key assumptions of the growth model include:
The natural mortality and recapture probability were estimated using a Bayesian individual state-space survival model (Thorley and Andrusak 2017) with three month intervals. The survival model incorporated natural and handling mortality, acoustic detections, inter-section movement, T-bar tag loss, spawning, recapture and reporting.
T-bar tag loss was estimated from the number of tags reported from recaught double-tagged individuals (Fabrizio et al. 1999).
Fork length at spawning was calculated from the measured length at
capture (\(L_{i,fi}\)) plus the length increment expected under a Von
Bertalanffy Growth Curve (Walters & Martell, 2004)
\[L_{i,t} = L_{i,fi} + (L_{∞} − L_{i,fi} )(1 − e^{(−0.25 \cdot k(t −fi))})\]
where the 0.25 in the exponent adjusts for the fact that there are four
time periods per year. Parameters k and LInf were estimated from the
growth model.
A fish was considered to have spawned in a given year if one of the following conditions were met: detected within a spawning river during the spawning period with a detection hiatus of at least 7 days; or detected within a spawning gate during the spawning period with a detection hiatus of at least 14 days. Cutoff values for minimum hiatus periods were determined by inspecting the distribution of in-lake hiatus periods for fish detected in spawning rivers and the distribution of all hiatus periods for suspected spawners within the spawning period.
In addition to assumptions 1 to 2 and 4 to 10, in Thorley and Andrusak (2017), the model also assumes that:
The optimal recapture rate (to maximize number of individuals caught) was calculated using a yield-per-recruit approach (Bison, O’Brien, and Martell 2003). Key assumptions include:
The optimal yield was also calculated with no slot limit.
.model {
bLInf ~ dnorm(900, 200^-2) T(0, )
bK ~ dnorm(0.2, 1^-2) T(0, )
bT0 ~ dnorm(0, 2^-1)
sLength ~ dnorm(0, 100^-2) T(0, )
for (i in 1:length(Length)) {
eLength[i] <- bLInf * (1 - exp(-bK * (Age[i] - bT0)))
Length[i] ~ dnorm(eLength[i], sLength^-2) T(0, )
}Block 1. Growth model description.
data {
int nannual;
int nObs;
vector[nObs] length;
vector[nObs] weight;
int annual[nObs];
}
parameters {
real bWeight;
real bLength;
real<lower=0> sWeightAnnual;
vector[nannual] bWeightAnnual;
real<lower=0> sWeight;
}
model {
vector[nObs] eWeight;
bWeight ~ normal(-11, 2);
bLength ~ normal(3, 1);
sWeightAnnual ~ normal(0, 1);
for (i in 1:nannual) {
bWeightAnnual[i] ~ normal(0, sWeightAnnual);
}
sWeight ~ normal(0, 1);
for(i in 1:nObs) {
eWeight[i] = exp(bWeight + bLength * log(length[i]) + bWeightAnnual[annual[i]]);
weight[i] ~ lognormal(log(eWeight[i]), sWeight);
}
}Block 2. Condition model description.
.model{
bMortality ~ dnorm(-3, 3^-2)
bMortalitySpawning ~ dnorm(0, 2^-2)
bSpawning ~ dnorm(0, 2^-2)
bSpawningLength ~ dnorm(0, 2^-2)
bTagLoss ~ dbeta(101, 1001)
bDetectedAlive ~ dbeta(1, 1)
bDieNear ~ dbeta(1, 1)
bDetectedDeadNear ~ dbeta(10, 1)
bDetectedDeadFar ~ dbeta(1, 10)
bMoved ~ dbeta(1, 1)
bRecaptured ~ dbeta(1, 1)
bReported ~ dbeta(10, 1)
sSpawningAnnual ~ dexp(1)
for (i in 1:nAnnual) {
bSpawningAnnual[i] ~ dnorm(0, sSpawningAnnual^-2)
}
bMortalitySeason[1] <- 0
for (i in 2:nSeason) {
bMortalitySeason[i] ~ dnorm(0, 2^-2)
}
sMortalityAnnual ~ dexp(1)
for (i in 1:nAnnual) {
bMortalityAnnual[i] ~ dnorm(0, sMortalityAnnual^-2)
}
for (i in 1:nCapture){
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 + bMortalitySeason[Season[i,PeriodCapture[i]]] + bMortalityAnnual[Annual[i,PeriodCapture[i]]]
eTagLoss[i,PeriodCapture[i]] <- bTagLoss
eMoved[i,PeriodCapture[i]] <- bMoved
eRecaptured[i,PeriodCapture[i]] <- bRecaptured
eReported[i,PeriodCapture[i]] <- bReported
InLake[i,PeriodCapture[i]] <- 1
Alive[i,PeriodCapture[i]] ~ dbern(1-eMortality[i,PeriodCapture[i]])
TBarTag100[i,PeriodCapture[i]] ~ dbern(1-eTagLoss[i,PeriodCapture[i]])
TBarTag10[i,PeriodCapture[i]] ~ dbern(1-eTagLoss[i,PeriodCapture[i]])
eDetected[i,PeriodCapture[i]] <- Alive[i,PeriodCapture[i]] * eDetectedAlive[i] + (1-Alive[i,PeriodCapture[i]]) * eDetectedDead[i]
Detected[i,PeriodCapture[i]] ~ dbern(Monitored[i,PeriodCapture[i]] * eDetected[i,PeriodCapture[i]])
Moved[i,PeriodCapture[i]] ~ dbern(Alive[i,PeriodCapture[i]] * Detected[i,PeriodCapture[i]] * eMoved[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]] * step(TBarTag100[i,PeriodCapture[i]] + TBarTag10[i,PeriodCapture[i]] - 1))
for(j in (PeriodCapture[i]+1):nPeriod) {
logit(eSpawning[i,j]) <- bSpawning + bSpawningLength * Length[i,j] + bSpawningAnnual[Annual[i,j]]
Spawned[i,j] ~ dbern(Alive[i,j-1] * Monitored[i,j] * SpawningPeriod[i,j] * eSpawning[i,j])
logit(eMortality[i,j]) <- bMortality + bMortalitySpawning * Spawned[i,j] + bMortalitySeason[Season[i,j]] + bMortalityAnnual[Annual[i,j]]
eDetected[i,j] <- Alive[i,j] * eDetectedAlive[i] + (1-Alive[i,j]) * eDetectedDead[i]
eTagLoss[i,j] <- bTagLoss
eMoved[i,j] <- bMoved
eRecaptured[i,j] <- bRecaptured
eReported[i,j] <- bReported
InLake[i,j] ~ dbern(InLake[i,j-1] * (1 - Recaptured[i,j-1]))
Alive[i,j] ~ dbern(Alive[i,j-1] * (1 - Recaptured[i,j-1]) * (1-eMortality[i,j]))
TBarTag100[i,j] ~ dbern(TBarTag100[i,j-1] * (1-Recaptured[i,j-1]) * (1-eTagLoss[i,j]))
TBarTag10[i,j] ~ dbern(TBarTag10[i,j-1] * (1-Recaptured[i,j-1]) * (1-eTagLoss[i,j]))
Detected[i,j] ~ dbern(InLake[i,j] * Monitored[i,j] * eDetected[i,j])
Moved[i,j] ~ dbern(Alive[i,j] * Detected[i,j] * eMoved[i,j])
Recaptured[i,j] ~ dbern(Alive[i,j] * eRecaptured[i,j])
Reported[i,j] ~ dbern(Recaptured[i,j] * eReported[i,j] * step(TBarTag100[i,j] + TBarTag10[i,j] - 1))
}
}Block 3. Survival model description.
The suggested citation for this analytic appendix is:
Thorley, J.L., Andrusak, G.F. and Amies-Galonski, E. (2022) Kaslo Bull Trout Productivity 2021. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/433821282.
The Kaslo River and Keen Creek, which is a tributary of the Kaslo River, are important Bull Trout spawning and rearing tributaries on Kootenay Lake. From 2012 to 2021, field crews have night-snorkeled these systems in the fall and recorded all Bull Trout less than 350 mm in length. Keen Creek was not snorkelled in 2020 or 2021 due to visibility issues. Snorkel and electrofishing marking crews have also captured and tagged juvenile Bull Trout for the snorkel crews to resight. Redd counts have been conducted in both systems since 2006, with the exception of 2020 and 2021, when Keen Creek was not surveyed. The primary goal of the current analyses is to answer the following questions:
What is the observer efficiency when night-snorkeling for juvenile Bull Trout in the Kaslo River and Keen Creek?
What are the numbers of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the stock (number of redds and eggs) and the resultant numbers of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 4.2.1 (R Core Team 2015).
Model parameters were estimated using Bayesian methods. The Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on Bayesian estimation the reader is referred to McElreath (2016).
Unless indicated otherwise, the Bayesian analyses used weakly informative normal or truncated half normal prior distributions (Kery and Schaub 2011, 36). 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 \(\hat{R} \leq\) 1.05 (Kery and Schaub 2011, 40) and \(\textrm{ESS} \geq 150\) for each of the monitored parameters (Kery and Schaub 2011, 61). Where \(\hat{R}\) is the potential scale reduction factor and \(\textrm{ESS}\) is the effective sample size.
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 can be considered a test of directionality. More specifically it indicates how surprising (in bits) it would be to discover that the true value of the parameter is in the opposite direction to the estimate. An s-value of 4.3 bits, which is equivalent to a p-value (Kery and Schaub 2011; Greenland and Poole 2013) of 0.05, indicates that the surprise would be equivalent to throwing 4.3 heads in a row.
Where relevant, model adequacy was confirmed by examination of residual plots for the full model(s).
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). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford et al. 2005).
The analyses were implemented using R version 4.2.1
(R Core Team 2015) and the
mbr family of packages.
The annual bias (inaccuracy) and error (imprecision) in observer’s fish length estimates when spotlighting (standing) and snorkeling were quantified from the divergence of their length distribution from the length distribution for JLT and SH (the two most experience snorkelers) in that year. More specifically, the length correction that minimised the Jensen-Shannon divergence (Lin 1991) provided a measure of the inaccuracy while the minimum divergence (the Jensen-Shannon divergence was calculated with log to base 2 which means it lies between 0 and 1) provided a measure of the imprecision.
All resighted fish with a tag were allocated to the closest unallocated marked fish (with the same colour tag) by fork length and distance. The marked fish were analysed using a Bayesian logistic regression model. The key assumption of the logistic regression model is that:
The preliminary analysis for observer efficiency indicated that system, observer, gradient, sinuosity, and river kilometre were not informative predictors.
Both systems were broken into 100 m sites by bank. The lineal density at each site was estimated using an over-dispersed Poisson Generalized Linear Mixed Model.
Key assumptions of the Bayesian GLMM include:
The preliminary analysis for density indicated that site sinuosity, gradient, and river kilometre were not informative predictors.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
The \(E_{K/2}\) Limit Reference Point (Mace 1994) (\(E_{0.5 R_{max}}\)) was calculated, corresponding to the stock (number of redds per kilometer) that produce 50% of the maximum recruitment (\(K\)).
The condition of adults was estimated using an allometric weight-length relationship. Key assumptions of the condition model include:
Female spawner fecundity was estimated using an allometric egg-weight relationship. The key assumption of the fecundity model is that:
The average length of spawners in each system and year was estimated using a linear regression. Key assumptions of the spawner length model include:
The total egg deposition in each year was estimated by
The stock-recruitment relationship between the total number eggs and the abundance of age-1 individuals was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
The \(E_{K/2}\) Limit Reference Point (Mace 1994) (\(E_{0.5 R_{max}}\)) was calculated, corresponding to the stock (number of eggs per 100 meters) that produce 50% of the maximum recruitment (\(K\)).
model{
bIntercept ~ dnorm(0, 2^-2)
bLength ~ dnorm(0, 2^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Block 1. Final model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, 1)
b0 ~ dnorm(0, 5^-2)
bRkm ~ dnorm(0, 2^-2)
bSystem[1] <- 0
for(i in 2:nSystem) {
bSystem[i] ~ dnorm(0, 2^2)
}
sSystemYear ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, sSystemYear^-2)
}
}
sSite ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
sDispersion ~ dnorm(0, 2^-2) T(0,)
for(i in 1:length(Count)) {
log(eDensity[i]) <- b0 + bSystem[System[i]] + bRkm * Rkm[i] + bSite[Site[i]] + bSystemYear[System[i],Year[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Block 2. The final model.
model {
log_bAlpha ~ dnorm(5, 5^-2)
log(bAlpha) <- log_bAlpha
for(i in 1:nSystem) {
log_bBeta[i] ~ dnorm(0, 5^-2)
log(bBeta[i]) <- log_bBeta[i]
}
log_sRecruits ~ dnorm(0, 5^-2)
log(sRecruits) <- log_sRecruits
for(i in 1:length(Recruits)){
eRecruits[i] <- bAlpha * Stock[i] / (1 + bBeta[System[i]] * Stock[i])
Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
}
bCarryingCapacity <- bAlpha / bBetaBlock 3. Final model.
model{
b0 ~ dnorm(6, 2^-2)
bLength ~ dnorm(3, 1^-2)
sWeight ~ dnorm(0, 1^-2) T(0,)
sYear ~ dnorm(0, 1^-2) T(0,)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nObs) {
log(eWeight[i]) <- b0 + (log(Length[i]) - log(500)) * bLength + bYear[Year[i]]
Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
}Block 4. Model description.
model{
b0 ~ dnorm(0, 2^-2)
bWeight ~ dnorm(1, 2^-2)
sFecundity ~ dnorm(0, 2^-2) T(0,)
for (i in 1:nObs) {
log(eFecundity[i]) <- b0 + (log(Weight[i]) - log(2300)) * bWeight
Fecundity[i] ~ dlnorm(log(eFecundity[i]), sFecundity^-2)
}Block 5. Model description.
model{
bLength ~ dnorm(650, 100^-2)
sLength ~ dnorm(0, 100^-2) T(0,)
sYear ~ dnorm(0, 100^-2) T(0,)
sSystem ~ dnorm(0, 100^-2) T(0,)
sYearSystem ~ dnorm(0, 100^-2) T(0,)
bSex ~ dbeta(1, 1)
bFemale ~ dnorm(0, 100^-2)
bBias ~ dnorm(0, 100^-2)
for (i in 1:nYear) {
bYear[i] ~ dnorm(0, sYear^-2)
}
for (i in 1:nSystem) {
bSystem[i] ~ dnorm(0, sSystem^-2)
}
for (i in 1:nYear) {
for (j in 1:nSystem) {
bYearSystem[i, j] ~ dnorm(0, sYearSystem^-2)
}
}
bCatchType[1] <- 0
for(i in 2:nCatchType) {
bCatchType[i] ~ dnorm(0, 100^-2)
}
for (i in 1:nObs) {
Female[i] ~ dbern(bSex)
eLength[i] <- bLength + bSystem[System[i]] + bFemale * Female[i] + bCatchType[CatchType[i]] + bYear[Year[i]] + bYearSystem[Year[i], System[i]] + bBias * Bias[i]
Length[i] ~ dnorm(eLength[i], sLength^-2)
}Block 6. Model description.
model {
bAlpha ~ dbeta(4, 49)
bK ~ dnorm(3, 1^-2)
bKPopulation ~ dnorm(0, 1^-2)
sRecruits ~ dexp(1)
for(i in 1:nObs){
log(eK[i]) <- bK + (Kaslo[i] * bKPopulation - (1-Kaslo[i]) * bKPopulation)
eBeta[i] <- bAlpha/eK[i]
eRecruits[i] <- bAlpha * Stock[i] / (1 + eBeta[i] * Stock[i])
Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
}Block 7. Final model.
The suggested citation for this analytic report is:
Thorley, J.L. & Hussein, N. (2022) Lower Columbia River Fish Population Indexing 2021. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1314086833.
In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.
The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:
What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River? What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River? The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.
The suggested citation for this analytic report is:
Thorley, J.L. & Hussein, N. (2021) Lower Columbia River Fish Population Indexing 2020. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/392669554.
In the mid 1990s BC Hydro began operating Hugh L. Keenleyside (HLK) Dam to reduce dewatering of Mountain Whitefish and Rainbow Trout eggs.
The primary goal of the Lower Columbia River Fish Population Indexing program is to answer two key management questions:
What are the abundance, growth rate, survival rate, body condition, age distribution, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
What is the effect of inter-annual variability in the Whitefish and Rainbow Trout flow regimes on the abundance, growth rate, survival rate, body condition, and spatial distribution of subadult and adult Whitefish, Rainbow Trout, and Walleye in the Lower Columbia River?
The inter-annual variability in the Whitefish and Rainbow Trout flow regimes was quantified in terms of the percent egg dewatering as greater flow variability is associated with more egg stranding.
The suggested citation for this analytic appendix is:
Thorley, J.L. & Amies-Galonski, E. (2021) Mica Dam Expansion Water Temperature and Fish Indexing Study 2020. A Poisson Consulting Analysis Appendix. URL: https://www.poissonconsulting.ca/f/955630057.
The Mica Tailrace Fish Indexing Study is a multi-year program to estimate the effects of the addition of two new turbines (Mica 5 and 6) on the ichyofauna and thermal regime in the 2.5 km of the Columbia River downstream of Mica Dam. A single year of fish indexing data (2008) was also available from a previous program. As per the Terms of Reference (TOR) the relative abundance, condition and spatial distribution of the fish populations was assessed.
Mica 5 became operational on January 28th 2015 and Mica 6 became operational on December 22nd 2015.
The fish and downstream temperature data were provided by the Ktunaxa Nation. The discharge and elevation data were queried from the Columbia Basin Hydrological Database.
The data were cleaned and tidied using R version 4.0.5 (R Core Team 2020).
Individuals were classified as fry (age-0), juvenile (age-1 and older subadults) or adult (sexually mature) based on the length cut-offs in Table 1.
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 (2016).
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 can be considered a test of directionality. More specifically it indicates how surprising (in bits) it would be to discover that the true value of the parameter is in the opposite direction to the estimate. An s-value of 4.3 bits, which is equivalent to a p-value (Kery and Schaub 2011; Greenland and Poole 2013) of 0.05, indicates that the surprise would be equivalent to throwing 4.3 heads in a row. The requirement that non-essential explanatory variables have s-values \(\geq\) 4.3 bits provides a useful model selection heuristic (Kery and Schaub 2011).
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). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent or n-fold change in the response variable) with 95% credible intervals (CIs, Bradford, Korman, and Higgins 2005).
The analyses were implemented using R version 4.0.5
(R Core Team 2020) and the
mbr family of packages.
The annual variation in condition (body weight when accounting for body length) was estimated from the boat and backpack electrofishing captures using a mass-length model (He et al. 2008).
Key assumptions of the condition model include:
Preliminary analyses indicated that site and day of the year were not informative predictors of condition.
The annual variation in relative abundance was estimated from the boat count and catch data using an over-dispersed Poisson model (Kery and Schaub 2011). Lineal densities are by kilometre of river (as opposed to kilometre of bank).
Key assumptions of the relative abundance model include:
Preliminary analyses indicated that site and discharge were not informative predictors of the lineal count (or catch) density.
The model does not distinguish between the abundance and observer efficiency, i.e., it estimates the count which is the product of the two. As such it is necessary to assume that changes in observer efficiency by year are negligible in order to interpret the estimates as relative abundance.
Climatic variation can cause large differences in annual temperatures. Consequently, we explored the data for an effect of the additional turbines on the difference in the water temperature between the right versus left bank and when moving downstream. All apparent systematic differences were within the accuracy of the temperature loggers (\(\pm 0.2^{\circ}\text{C}\)).
In 2020 the BC Hydro thermistor data for the Mica Dam forebay were analysed using Maximum Likelihood (Millar 2011) and TMB (Kristensen et al. 2016) to determine the extent to which discharge from the turbines and the air temperature influences the water temperature in the tailrace.
Key assumptions of the hourly temperature model include:
Despite the incorporation of a first order moving average process the residual variation in the hourly tailrace water temperature was strongly autocorrelated. Consequently CIs are not provided for the coefficients or estimates.
.model{
bAlpha ~ dnorm(5, 2^-2)
bBeta ~ dnorm(3, 2^-2)
bAlphaPeriod[1] <- 0
for(i in 2:nPeriod) {
bAlphaPeriod[i] ~ dnorm(0, 2^-2)
}
sAlphaAnnual ~ dnorm(0, 2^-2) T(0,)
for(i in 1:nAnnual) {
bAlphaAnnual[i] ~ dnorm(0, sAlphaAnnual^-2)
}
sWeight ~ dnorm(0, 2^-2) T(0,)
for (i in 1:length(Weight)) {
eAlpha[i] <- bAlpha + bAlphaPeriod[Period[i]] + bAlphaAnnual[Annual[i]]
eBeta[i] <- bBeta
log(eWeight[i]) <- eAlpha[i] + eBeta[i] * Length[i]
Weight[i] ~ dlnorm(log(eWeight[i]), sWeight^-2)
}Block 1.
.model{
bEfficiencyVisitType[1] <- 0
for (i in 2:nVisitType) {
bEfficiencyVisitType[i] ~ dnorm(0, 1^-2)
}
bDensity ~ dnorm(0, 5^-2)
bDensityPeriod[1] <- 0
for(i in 2:nPeriod) {
bDensityPeriod[i] ~ dnorm(0, 2^-2)
}
sDensityAnnual ~ dnorm(0, 2^-2) T(0, )
for(i in 1:nAnnual) {
bDensityAnnual[i] ~ dnorm(0, sDensityAnnual^-2)
}
sDispersion ~ dnorm(0, 2^-2) T(0, )
for (i in 1:length(Annual)) {
log(eEfficiency[i]) <- bEfficiencyVisitType[VisitType[i]]
log(eDensity[i]) <- bDensity + bDensityPeriod[Period[i]] + bDensityAnnual[Annual[i]]
eAbundance[i] <- eDensity[i] * SiteLength[i] / 2
eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)
Count[i] ~ dpois(eAbundance[i] * eEfficiency[i] * eDispersion[i])
}Block 2.
.#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Tailrace);
DATA_VECTOR(Intake);
DATA_VECTOR(Forebay);
DATA_VECTOR(Air);
DATA_VECTOR(Unit1);
DATA_VECTOR(Unit2);
DATA_VECTOR(Unit3);
DATA_VECTOR(Unit4);
DATA_VECTOR(Unit5);
DATA_VECTOR(Unit6);
DATA_VECTOR(Spill);
DATA_IVECTOR(Consecutive);
PARAMETER(sDepthUnit);
PARAMETER_VECTOR(bDepthUnit);
PARAMETER(bDepthSpill0);
PARAMETER_VECTOR(bDepthTurbineUnit);
PARAMETER(bRelease);
PARAMETER(bReleaseAir);
PARAMETER(bRho);
PARAMETER(sTailrace);
vector<Type> eTurbineUnit1 = Tailrace;
vector<Type> eTurbineUnit2 = Tailrace;
vector<Type> eTurbineUnit3 = Tailrace;
vector<Type> eTurbineUnit4 = Tailrace;
vector<Type> eTurbineUnit5 = Tailrace;
vector<Type> eTurbineUnit6 = Tailrace;
vector<Type> eDischarge = Tailrace;
vector<Type> eDepthTurbineUnit1 = Tailrace;
vector<Type> eDepthTurbineUnit2 = Tailrace;
vector<Type> eDepthTurbineUnit3 = Tailrace;
vector<Type> eDepthTurbineUnit4 = Tailrace;
vector<Type> eDepthTurbineUnit5 = Tailrace;
vector<Type> eDepthTurbineUnit6 = Tailrace;
vector<Type> eDepthSpill = Tailrace;
vector<Type> eStratification = Tailrace;
vector<Type> eEffectTurbineUnit1 = Tailrace;
vector<Type> eEffectTurbineUnit2 = Tailrace;
vector<Type> eEffectTurbineUnit3 = Tailrace;
vector<Type> eEffectTurbineUnit4 = Tailrace;
vector<Type> eEffectTurbineUnit5 = Tailrace;
vector<Type> eEffectTurbineUnit6 = Tailrace;
vector<Type> eEffectSpill = Tailrace;
vector<Type> eEffectDischarge = Tailrace;
vector<Type> eReleaseEffect = Tailrace;
vector<Type> eOutput = Tailrace;
vector<Type> eRelease = Tailrace;
vector<Type> eTailrace = Tailrace;
Type nll = 0.0;
eTailrace(0) = Tailrace(0);
for(int i = 0; i < 6; i++){
nll -= dnorm(bDepthUnit(i), Type(0), exp(sDepthUnit), true);
for(int i = 1; i < Tailrace.size(); i++){
eTurbineUnit1(i) = Unit1(i) + 0.0001;
eTurbineUnit2(i) = Unit2(i) + 0.0001;
eTurbineUnit3(i) = Unit3(i) + 0.0001;
eTurbineUnit4(i) = Unit4(i) + 0.0001;
eTurbineUnit5(i) = Unit5(i) + 0.0001;
eTurbineUnit6(i) = Unit6(i) + 0.0001;
eDischarge(i) = eTurbineUnit1(i) + eTurbineUnit2(i) + eTurbineUnit3(i) + eTurbineUnit4(i) + eTurbineUnit5(i) + eTurbineUnit6(i) + Spill(i);
eDepthTurbineUnit1(i) = bDepthUnit(0) + bDepthTurbineUnit(0) * eTurbineUnit1(i);
eDepthTurbineUnit2(i) = bDepthUnit(1) + bDepthTurbineUnit(1) * eTurbineUnit2(i);
eDepthTurbineUnit3(i) = bDepthUnit(2) + bDepthTurbineUnit(2) * eTurbineUnit3(i);
eDepthTurbineUnit4(i) = bDepthUnit(3) + bDepthTurbineUnit(3) * eTurbineUnit4(i);
eDepthTurbineUnit5(i) = bDepthUnit(4) + bDepthTurbineUnit(4) * eTurbineUnit5(i);
eDepthTurbineUnit6(i) = bDepthUnit(5) + bDepthTurbineUnit(5) * eTurbineUnit6(i);
eDepthSpill(i) = bDepthSpill0;
eStratification(i) = (Forebay(i) - Intake(i)) / (707 - 690) ;
eEffectTurbineUnit1(i) = eDepthTurbineUnit1(i) * eStratification(i);
eEffectTurbineUnit2(i) = eDepthTurbineUnit2(i) * eStratification(i);
eEffectTurbineUnit3(i) = eDepthTurbineUnit3(i) * eStratification(i);
eEffectTurbineUnit4(i) = eDepthTurbineUnit4(i) * eStratification(i);
eEffectTurbineUnit5(i) = eDepthTurbineUnit5(i) * eStratification(i);
eEffectTurbineUnit6(i) = eDepthTurbineUnit6(i) * eStratification(i);
eEffectSpill(i) = eDepthSpill(i) * eStratification(i);
eEffectDischarge(i) = (eEffectTurbineUnit1(i) * eTurbineUnit1(i) + eEffectTurbineUnit2(i) * eTurbineUnit2(i) + eEffectTurbineUnit3(i) * eTurbineUnit3(i) + eEffectTurbineUnit4(i) * eTurbineUnit4(i) + eEffectTurbineUnit5(i) * eTurbineUnit5(i) + eEffectTurbineUnit6(i) * eTurbineUnit6(i) + eEffectSpill(i) * Spill(i)) / eDischarge(i);
eReleaseEffect(i) = exp(bRelease + bReleaseAir * (Air(i) - Intake(i)));
eOutput(i) = Intake(i) + eEffectDischarge(i);
eRelease(i) = eOutput(i) + eReleaseEffect(i);
eTailrace(i) = eRelease(i) + Consecutive(i) * invlogit(bRho) * (Tailrace(i-1) - eTailrace(i-1));
nll -= dnorm(Tailrace(i), eTailrace(i), exp(sTailrace), true);
return nll;Block 3. Model description.