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:
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:
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 appendix is:
Thorley, J.L., Andrusak, G.F., Amies-Galonski, E. and Dalgarno, S.I.J. (2020) Kaslo Bull Trout Productivity 2020. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/1283025452.
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 2020, field crews have night-snorkeled these systems in the fall and recorded all juvenile Bull Trout less than 350 mm in length. Keen Creek was not snorkelled in 2020 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, in which 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 densities of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 4.0.3 (R Core Team 2015).
Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML and Bayesian estimation the reader is referred to Millar (2011) and McElreath (2016), respectively.
Unless indicated otherwise, the Bayesian analyses used uninformative 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.1\) (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.
Model averaging and selection was implemented using an information theoretic approach (Burnham and Anderson 2002). ML models were compared using the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The support for fixed terms in models with identical random (Kery and Schaub 2011, 75) effects was assessed using the Akaike’s weights (\(w_i\)) (Burnham and Anderson 2002). The best supported model was then refitted as its Bayesian equivalent.
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.0.3
(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 logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:
Preliminary analysis indicated that river kilometre received similar support to watershed area.
The final Bayesian logistic regression assumed that:
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 full Maximum Likelihood GLMM include:
Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area. Preliminary analysis also indicated that autocorrelation (modeled as an AR1 process) was not significant.
The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
.#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(Observer);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bObserver);
vector<Type> eObserved = Observed;
int n = Observed.size();
Type nll = 0.0;
for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bObserver(Observer(i)) + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;Block 1. Full model.
model{
bIntercept ~ dnorm(0, 5^-2)
bLength ~ dnorm(0, 5^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Block 2. Final model.
.#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
DATA_FACTOR(System);
DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);
PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);
Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);
ADREPORT(sSite);
ADREPORT(sDispersion);
vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));
Type nll = 0.0;
for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);
for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);
eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);
nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);
return nll;Block 3. The full model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, 5^-2)
}
}
log_sSite ~ dnorm(0, 5^-2)
log(sSite) <- log_sSite
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
bRkm ~ dnorm(0, 5^-2)
log_sDispersion ~ dnorm(0, 5^-2)
log(sDispersion) <- log_sDispersion
for(i in 1:length(Count)) {
log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bRkm * Rkm[i] + bSite[Site[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Block 4. 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 5. Final model.
The suggested citation for this analytic appendix is:
Thorley, J.L., Andrusak, G.F., Amies-Galonski, E. and Dalgarno, S.I.J. (2019) Kaslo Bull Trout Productivity 2019. A Poisson Consulting Analytic Appendix. URL: https://www.poissonconsulting.ca/f/168922378.
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 2019, field crews have night-snorkeled both systems in the fall and recorded all juvenile Bull Trout less than 350 mm in length. 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.
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 densities of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 3.6.1 (R Core Team 2015).
Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML and Bayesian estimation the reader is referred to Millar (2011) and McElreath (2016), respectively.
Unless indicated otherwise, the Bayesian analyses used uninformative 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.1\) (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, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). For ML models, the point estimate is the MLE, the standard deviation is the standard error, the z-score is \(\mathrm{sd}/\mathrm{MLE}\) and the 95% CLs are the \(\mathrm{MLE} \pm 1.96 \cdot \mathrm{sd}\). For Bayesian models, the estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{sd}/\mathrm{mean}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.
Model averaging and selection was implemented using an information theoretic approach (Burnham and Anderson 2002). ML models were compared using the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The support for fixed terms in models with identical random (Kery and Schaub 2011, 75) effects was assessed using the Akaike’s weights (\(w_i\)) (Burnham and Anderson 2002). The best supported model was then refitted as its Bayesian equivalent.
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 3.6.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 logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:
Preliminary analysis indicated that river kilometre received similar support to watershed area.
The final Bayesian logistic regression assumed that:
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 full Maximum Likelihood GLMM include:
Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area. Preliminary analysis also indicated that autocorrelation (modeled as an AR1 process) was not significant.
The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
.#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(Observer);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bObserver);
vector<Type> eObserved = Observed;
int n = Observed.size();
Type nll = 0.0;
for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bObserver(Observer(i)) + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;Block 1. Full model.
model{
bIntercept ~ dnorm(0, 5^-2)
bLength ~ dnorm(0, 5^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Block 2. Final model.
.#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
DATA_FACTOR(System);
DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);
PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);
Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);
ADREPORT(sSite);
ADREPORT(sDispersion);
vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));
Type nll = 0.0;
for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);
for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);
eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);
nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);
return nll;Block 3. The full model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, 5^-2)
}
}
log_sSite ~ dnorm(0, 5^-2)
log(sSite) <- log_sSite
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
bRkm ~ dnorm(0, 5^-2)
log_sDispersion ~ dnorm(0, 5^-2)
log(sDispersion) <- log_sDispersion
for(i in 1:length(Count)) {
log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bRkm * Rkm[i] + bSite[Site[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Block 4. 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 5. Final model.
The suggested citation for this analytic report is:
Thorley, J.L., Andrusak, G.F., Dalgarno, S.I.J., Amies-Galonski, E. (2018) Kaslo Bull Trout Productivity 2018. A Poisson Consulting Analysis Report. URL: https://www.poissonconsulting.ca/f/1920399815.
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 2018, field crews have night-snorkeled both systems in the fall and recorded all juvenile Bull Trout less than 350 mm in length. 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.
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 densities of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 3.5.1 (R Core Team 2015).
Table 1. Bull Trout length cutoffs by age-class.
| Age Class | Length (mm) |
|---|---|
| Age-0 | 30 - 90 |
| Age-1 | 91 - 140 |
| Age-2+ | 141 - 350 |
Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML and Bayesian estimation the reader is referred to Millar (2011) and McElreath (2016), respectively.
Unless indicated otherwise, the Bayesian analyses used uninformative 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.1\) (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, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). For ML models, the point estimate is the MLE, the standard deviation is the standard error, the z-score is \(\mathrm{sd}/\mathrm{MLE}\) and the 95% CLs are the \(\mathrm{MLE} \pm 1.96 \cdot \mathrm{sd}\). For Bayesian models, the estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{mean}/\mathrm{sd}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.
Model averaging and selection was implemented using an information theoretic approach (Burnham and Anderson 2002). ML models were compared using the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The support for fixed terms in models with identical random (Kery and Schaub 2011, 75) effects was assessed using the Akaike’s weights (\(w_i\)) (Burnham and Anderson 2002). The best supported model was then refitted as its Bayesian equivalent.
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 3.5.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 logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:
Preliminary analysis indicated that river kilometre received similar support to watershed area.
The final Bayesian logistic regression assumed that:
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 full Maximum Likelihood GLMM include:
Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area. Preliminary analysis also indicated that autocorrelation (modeled as an AR1 process) was not significant.
The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(Observer);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bObserver);
vector<Type> eObserved = Observed;
int n = Observed.size();
Type nll = 0.0;
for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bObserver(Observer(i)) + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;Template 1. Full model.
model{
bIntercept ~ dnorm(0, 5^-2)
bLength ~ dnorm(0, 5^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Template 2. Final model.
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
DATA_FACTOR(System);
DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);
PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);
Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);
ADREPORT(sSite);
ADREPORT(sDispersion);
vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));
Type nll = 0.0;
for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);
for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);
eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);
nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);
return nll;Template 3. The full model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, 5^-2)
}
}
log_sSite ~ dnorm(0, 5^-2)
log(sSite) <- log_sSite
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
bRkm ~ dnorm(0, 5^-2)
log_sDispersion ~ dnorm(0, 5^-2)
log(sDispersion) <- log_sDispersion
for(i in 1:length(Count)) {
log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bRkm * Rkm[i] + bSite[Site[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Template 4. 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 / bBetaTemplate 5. Final model.
The suggested citation for this analytic report is:
Thorley, J.L., Andrusak, G.F. and Dalgarno, S.I.J. (2018) Kaslo Bull Trout Productivity 2017. A Poisson Consulting Analysis Report. URL: https://www.poissonconsulting.ca/f/443864882.
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 2017, field crews have night-snorkeled both systems in the fall and recorded all juvenile Bull Trout between 30 mm and 350 mm in length. 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.
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 densities of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?
The data were cleaned, tidied and databased using R version 3.4.4 (R Core Team 2015).
Table 1. Bull Trout length cutoffs by age-class.
| Age Class | Length (mm) |
|---|---|
| Age-0 | 30 - 90 |
| Age-1 | 91 - 150 |
| Age-2+ | 151 - 350 |
Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML and Bayesian estimation the reader is referred to Millar (2011) and McElreath (2016), respectively.
Unless indicated otherwise, the Bayesian analyses used uninformative 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.1\) (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, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). For ML models, the point estimate is the MLE, the standard deviation is the standard error, the z-score is \(\mathrm{sd}/\mathrm{MLE}\) and the 95% CLs are the \(\mathrm{MLE} \pm 1.96 \cdot \mathrm{sd}\). For Bayesian models, the estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{mean}/\mathrm{sd}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.
Model averaging and selection was implemented using an information theoretic approach (Burnham and Anderson 2002). ML models were compared using the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The support for fixed terms in models with identical random (Kery and Schaub 2011, 75) effects was assessed using the Akaike’s weights (\(w_i\)) which were also used for model averaging (Burnham and Anderson 2002).
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 3.4.4
(R Core Team 2015) and the
mbr family of packages.
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 logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:
Preliminary analysis indicated that river kilometre received similar support to watershed area.
The final Bayesian logistic regression assumed that:
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 full Maximum Likelihood GLMM include:
Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area.
The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
PARAMETER(bIntercept);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
vector<Type> eObserved = Observed;
int n = Observed.size();
Type nll = 0.0;
for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bIntercept + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;Template 1. Full model.
model{
bIntercept ~ dnorm(0, 5^-2)
bLength ~ dnorm(0, 5^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Template 2. Final model.
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);
DATA_FACTOR(System);
DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);
PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);
PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);
Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);
ADREPORT(sSite);
ADREPORT(sDispersion);
vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));
Type nll = 0.0;
for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);
for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);
eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);
nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);
return nll;Template 3. The full model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, 5^-2)
}
}
log_sSite ~ dnorm(0, 5^-2)
log(sSite) <- log_sSite
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
bRkm ~ dnorm(0, 5^-2)
log_sDispersion ~ dnorm(0, 5^-2)
log(sDispersion) <- log_sDispersion
for(i in 1:length(Count)) {
log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bRkm * Rkm[i] + bSite[Site[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Template 4. 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 / bBetaTemplate 5. Final model.
The suggested citation for this analytic report is:
Thorley, J.L., Andrusak, G.F. and Hogan P.M. (2017) Kaslo Bull Trout Productivity Analysis 2016. A Poisson Consulting Analysis Report. URL: https://www.poissonconsulting.ca/f/1827953676.
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 2016, field crews have night-snorkeled both systems in the fall and recorded all juvenile Bull Trout between 30 mm and 350 mm in length. 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.
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 densities of age-1 Bull Trout in the Kaslo River and Keen Creek?
What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?
The data were provided in the form of clean and tidy Excel spreadsheets by Gary Pavan and prepared for analysis using R version 3.3.2 (R Core Team 2015).
Table 1. Bull Trout length cutoffs by age-class.
| Age Class | Length (mm) |
|---|---|
| Age-0 | 30 - 90 |
| Age-1 | 91 - 150 |
| Age-2+ | 151 - 350 |
Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML estimation the reader is referred to Millar (2011). For additional information on Bayesian modelling in the BUGS language, of which JAGS uses a dialect, the reader is referred to Kery and Schaub (2011).
Unless indicated otherwise, the Bayesian analyses used uninformative normal prior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from 2,000 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of four chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that \(\hat{R} < 1.1\) (Kery and Schaub 2011, 40) for each of the monitored parameters (Kery and Schaub 2011, 61).
The parameters are summarised in terms of the point estimate, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). For ML models, the point estimate is the MLE, the standard deviation is the standard error, the z-score is \(\mathrm{sd}/\mathrm{MLE}\) and the 95% CLs are the \(\mathrm{MLE} \pm 1.96 \cdot \mathrm{sd}\). A p-value of 0.05 indicates that the lower or upper 95% CL is 0. For Bayesian models, the estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{mean}/\mathrm{sd}\) and the 95% CLs are the 2.5th and 97.5th percentiles.
The support for fixed terms in ML models with identical random (Kery and Schaub 2011, 75) effects was assessed using Akaike’s weights (\(w_i\)) calculated from the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The Akaike’s weights were used for model averaging (Burnham and Anderson 2002). The fixed terms with low support were dropped from the final Bayesian model.
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 3.3.2
(R Core Team 2015) and the tmbr and jmbr packages.
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 logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:
The final Bayesian logistic regression assumed that:
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 full Maximum Likelihood GLMM include:
The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.
The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Area);
PARAMETER(bIntercept);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bArea);
vector<Type> eObserved = Observed;
int n = Observed.size();
Type nll = 0.0;
for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bIntercept + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bArea * Area(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;Template 1. Full model.
model{
bIntercept ~ dnorm(0, 5^-2)
bLength ~ dnorm(0, 5^-2)
for(i in 1:length(Observed)){
logit(eObserved[i]) <- bIntercept + bLength * Length[i]
Observed[i] ~ dbern(eObserved[i])
}Template 2. Final model.
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Area);
DATA_FACTOR(System);
DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);
PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bArea);
PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);
Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);
ADREPORT(sSite);
ADREPORT(sDispersion);
vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));
Type nll = 0.0;
for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);
}
for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);
eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bArea * Area(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);
nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);
}
return nll;
}Template 3. The full model.
model{
bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )
for(i in 1:nSystem) {
for(j in 1:nYear) {
bSystemYear[i,j] ~ dnorm(0, 5^-2)
}
}
log_sSite ~ dnorm(0, 5^-2)
log(sSite) <- log_sSite
for(i in 1:nSite) {
bSite[i] ~ dnorm(0, sSite^-2)
}
bArea ~ dnorm(0, 5^-2)
log_sDispersion ~ dnorm(0, 5^-2)
log(sDispersion) <- log_sDispersion
for(i in 1:length(Count)) {
log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bArea * Area[i] + bSite[Site[i]]
eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
}Template 4. 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 / bBetaTemplate 5. Final model.