Abundance

Data Preparation

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.

Statistical Analysis

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.

Model Descriptions

Model Templates

Observer Efficiency

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

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])
  }

Lineal Density

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

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])
  }

Stock-Recruiment

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 / bBeta

Data Preparation

The data were provided in the form of an Access database. The data were rearchived in an SQLite database and prepared for analysis using R version 3.5.0 (R Core Team 2018).

Statistical Analysis

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 normal and uniform prior distributions that did not constrain the posteriors (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 (Brooks et al. 2011).

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). 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.

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, Korman, and Higgins 2005).

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

Model Descriptions

Density

The lineal density and observer efficiency for Rainbow \(\geq\) 30 cm in length were estimated from the snorkel count and mark-resighting data for the Kettle River using an overdispersed Poisson mark-recapture model.

Key assumptions of the density model include:

• Density varies randomly by year and site.
• The number of resighted individuals is described by a binomial distribution where the probability of resighting is the observer efficiency.
• The expected number of observed individuals is the product of the density, site length, proportion of the site surveyed and the observer efficiency.
• The number of observed individuals is described by a gamma-Poisson distribution.

Model Templates

Density

model {
  bEfficiency ~ dunif(0, 1)

  bDensity ~ dnorm(5, 2^-2)

  sDensityAnnual ~ dunif(0, 1)
  for (yr in 1:nAnnual) {
    bDensityAnnual[yr] ~ dnorm(0, sDensityAnnual^-2)
  }

  sDensitySite ~ dunif(0, 2)
  for (i in 1:nSite) {
    bDensitySite[i] ~ dnorm(0, sDensitySite^-2)
  }

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

    eEfficiency[i] <- bEfficiency
    eResightingEfficiency[i] <- min(eEfficiency[i], Marked[i])

    eAdjustedMarked[i] <- max(Marked[i], 1)
    Resighted[i] ~ dbin(eResightingEfficiency[i] * ProportionSurveyed[i], eAdjustedMarked[i])

    log(eDensity[i]) <- bDensity + bDensityAnnual[Annual[i]] + bDensitySite[Site[i]]
    eAbundance[i] <- eDensity[i] * SiteLength[i]
    eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)
    Count[i] ~ dpois(eAbundance[i] * eEfficiency[i] * ProportionSurveyed[i] * eDispersion[i])
  }
..

Template 1. The model description.

Data Preparation

The age-1 Rainbow Trout abundance data were provided by the Fish Population Indexing Program. The remaining fish data were collected by Mountain Water Research and imported to an SQLite database.

The study area was divided into three sections: the LCR above the LKR, the LKR and the LCR below the LKR. Redd and spawner counts upstream of Norns Creek Fan and downstream of Genelle were excluded from the section totals because they constitute less than 0.1% of the total count and were not surveyed in all years. The redd and spawner counts for the Right Upstream Bank above Robson Bridge were also excluded as they appear to be primarily driven by viewing conditions (and constitute less than 2.5% of the total). Viewing conditions were classified as Good or Poor. If information on the viewing conditions was not available a decline in the redd count of more than one third of the cumulative maximum count for a particular section was assumed to be caused by poor viewing conditions.

The data were prepared for analysis using R version 3.4.4 (R Core Team 2017).

Statistical Analysis

Model parameters were estimated using Bayesian methods. The Bayesian estimates were produced using JAGS (Plummer 2015) and Stan (Carpenter et al. 2017). 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). For additional information on Bayesian modelling in the Stan language the reader is referred to Stan Development Team (2017).

Unless indicated otherwise, the Bayesian analyses used uninformative uniform or normal prior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from 1,500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of three chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that \(\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% credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). A p-value of 0.05 indicates that the lower or upper 95% CL is 0. 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 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, Korman, and Higgins 2005).

The analyses were implemented using R version 3.4.4 (R Core Team 2017) and the jmbr and smbr packages.

Model Descriptions

Stock-Recruitment

The relationship between the number of spawners and the resultant number of age-1 fish was estimated using a Beverton-Holt stock-recruitment model (Walters and Martell 2004):

\[ R = \frac{\alpha \cdot S}{1 + \beta \cdot S} \quad,\]

where \(S\) is the spawners (stock), \(R\) is the recruits, \(\alpha\) is the recruits per spawner at low density and \(\beta\) determines the density-dependence.

Key assumptions of the stock-recruitment model include:

• The recruits per spawner at low density (\(\alpha\)) is normally distributed with a mean of 90 and a SD of 50.
• The recruits per spawner varies with the percent of redds dewatered.
• The residual variation in the number of recruits is log-normally distributed.

The mean of 90 for \(\alpha\) was based on an average of 2,900 eggs per female spawner, a 50:50 sex ratio, 50% egg survival, 50% post-emergence fall survival, 50% overwintering survival and 50% summer survival.

The carrying capacity is \(\alpha / \beta\).

Model Templates

Area-Under-The-Curve

data {
  int<lower=0> nObs;
  int<lower=0> nSection;
  int<lower=0> nYear;

  int Year[nObs];
  int Section[nObs];
  real Doy[nObs];
  int Spawners[nObs];
  int Redds[nObs];
parameters {
  vector<lower=5,upper=10>[nSection] bSpawnerAbundanceSection;

  real<lower=0,upper=2> sSpawnerAbundanceYear;
  vector[nYear] bSpawnerAbundanceYear;
  real<lower=0,upper=2> sSpawnerAbundanceSectionYear;
  vector[nSection * nYear] bSpawnerAbundanceSectionYear;

  real<lower=0.8,upper=1.0> bSpawnerObserverEfficiency;
  real<lower=0.0,upper=2.0> bReddObserverEfficiency;

  real<lower=1, upper=2> bReddPerSpawner;

  real<lower=14, upper=21> bSpawnerResidenceTime;
  real<lower=30, upper=40> bReddResidenceTime;

  real<lower=100, upper=150> bPeakSpawnerArrivalTiming;

  real<lower=0,upper=21> sPeakSpawnerArrivalTimingYear;
  vector[nYear] bPeakSpawnerArrivalTimingYear;

  real<lower=log(10), upper=log(50)> bSpawnerArrivalDuration;

  real<lower=0,upper=2> sSpawnerArrivalDurationSectionYear;
  vector[nSection * nYear] bSpawnerArrivalDurationSectionYear;

  real<lower=0, upper=2> bDispersionRedds;
  real<lower=0, upper=2> bDispersionSpawners;
model {
  vector[nObs] eSpawnerAbundance;
  vector[nObs] ePeakSpawnerArrivalTiming;
  vector[nObs] eSpawnerArrivalDuration;

  vector[nObs] eRedds;
  vector[nObs] eSpawners;

  bSpawnerAbundanceYear ~ normal(0, sSpawnerAbundanceYear);
  bSpawnerAbundanceSectionYear ~ normal(0, sSpawnerAbundanceSectionYear);

  bPeakSpawnerArrivalTimingYear ~ normal(0, sPeakSpawnerArrivalTimingYear);
  bSpawnerArrivalDurationSectionYear ~ normal(0, sSpawnerArrivalDurationSectionYear);

  for (i in 1:nObs) {
    eSpawnerAbundance[i] = exp(bSpawnerAbundanceSection[Section[i]] + bSpawnerAbundanceYear[Year[i]] + bSpawnerAbundanceSectionYear[((Section[i] - 1) * Year[i]) + Year[i]]);

    ePeakSpawnerArrivalTiming[i] = bPeakSpawnerArrivalTiming + bPeakSpawnerArrivalTimingYear[Year[i]];

    eSpawnerArrivalDuration[i] = exp(bSpawnerArrivalDuration + bSpawnerArrivalDurationSectionYear[((Section[i] - 1) * Year[i]) + Year[i]]);

    eRedds[i] = eSpawnerAbundance[i] * bReddPerSpawner * bReddObserverEfficiency * (normal_cdf(Doy[i], ePeakSpawnerArrivalTiming[i], eSpawnerArrivalDuration[i]) - normal_cdf(Doy[i], ePeakSpawnerArrivalTiming[i] + bReddResidenceTime, eSpawnerArrivalDuration[i]));

    eSpawners[i] = eSpawnerAbundance[i] * bSpawnerObserverEfficiency * (normal_cdf(Doy[i], ePeakSpawnerArrivalTiming[i], eSpawnerArrivalDuration[i]) - normal_cdf(Doy[i], ePeakSpawnerArrivalTiming[i] + bSpawnerResidenceTime, eSpawnerArrivalDuration[i]));
  }
  Redds ~ neg_binomial_2(eRedds, 1/bDispersionRedds);
  Spawners ~ neg_binomial_2(eSpawners, 1/bDispersionSpawners);
..

Template 1. AUC model.

Stock-Recruitment

model {

  bAlpha ~ dnorm(90, 50^-2) T(1,)
  bAlphaDewatered ~ dnorm(0, 2^-2)
  bBeta ~ dnorm(-5, 5^-2)
  sRecruits ~ dunif(0, 2)

  for(i in 1:length(Stock)) {
    log(eAlpha[i]) <- log(bAlpha) + bAlphaDewatered * Dewatered[i]
    log(eBeta[i]) <- bBeta
    eRecruits[i] <- eAlpha[i] * Stock[i] / (1 + eBeta[i] * Stock[i])
    Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
  }
..

Template 2.

Data Preparation

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.

Statistical Analysis

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.

Model Descriptions

Model Templates

Observer Efficiency

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

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])
  }

Lineal Density

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

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])
  }

Stock-Recruiment

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 / bBeta

Data Preparation

The data were provided by BC Hydro, Ministry of Forests, Lands and Natural Resource Operations, Water Survey of Canada and Mountain Waters Research. The data were prepared for analysis using R version 3.4.4 (R Core Team 2017). The total spawner abundance at Gerrard was assumed to be the observed peak spawner count multiplied 3.08 (Irvine 1978).

Statistical Analysis

Model parameters were estimated using Bayesian methods (McElreath 2016).

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% credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). 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.

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, Korman, and Higgins 2005).

The analyses were implemented using R version 3.4.4 (R Core Team 2017), JAGS (Plummer 2015), and the mbr family of packages.

Model Descriptions

Model Templates

Redd Surveys

model{
    bAbundance ~ dnorm(0, 5^-2)
    bTiming ~ dunif(70, 140)
    bDuration ~ dunif(7, 28)

    sAbundanceYear ~ dunif(0, 5)
    sTimingYear ~ dunif(0, 40)
    sDurationYear ~ dunif(0, 20)
    for(i in 1:nYear){
      bTimingYear[i] ~ dnorm (0, sTimingYear^-2)
      bAbundanceYear[i] ~ dnorm (0, sAbundanceYear^-2)
      bDurationYear[i] ~ dnorm (0, sDurationYear^-2)
    }

    sRedds ~ dunif(0, 20)
    for(i in 1:nObs){
      log(eAbundance[i]) <- bAbundance + bAbundanceYear[Year[i]]
      eTiming[i] <- bTiming + bTimingYear[Year[i]]
      eDuration[i] <- bDuration + bDurationYear[Year[i]]

      eRedds[i] <- phi((Doy[i] - eTiming[i]) / eDuration[i]) * eAbundance[i]
      Redds[i] ~ dnorm(eRedds[i], sRedds^-2)
    }
  }

Template 1. The model description.

Data Preparation

The 2017 survey data were provided by LGL Limited in the form of Excel tables. They were imported into the existing SQLite database. The data were clean and tidied (Wickham 2014) using R version 3.4.3 (R Core Team 2016). When there were repeated counts for a channel the higher of the two values was choosen. The mean daily discharge and water temperature values at Duncan River below Lardeau (DRL) were extracted from BC Hydro’s environmental database for the Kootenays, which is maintained by Poisson Consulting. The annual abundances of Kokanee spawning in Meadow Creek and the Lardeau River were provided by the Ministry of Forests, Lands and Natural Resource Operations.

Statistical Analysis

Model parameters were estimated using Bayesian methods. The 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 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). 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.

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, Korman, and Higgins 2005).

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

Model Descriptions

Model Templates

Observer Efficiency

model{

  bEfficiency ~ dunif(0, 3)
  sDispersion ~ dunif(0, 5)

  for (i in 1:nObs){
    eAerial[i] <- Ground[i] * bEfficiency
    eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
    Aerial[i] ~ dpois(eAerial[i] * eDispersion[i])
  }
..

Template 1. Model definition.

Area-Under-The-Curve

model{
  bAbundance ~ dnorm(10, 5^-2)
  bPeakTiming ~ dunif(-14, 14)
  bDuration ~ dnorm(2, 1^-2)

  bEfficiency ~ dnorm(1.40, 0.38^-2) T(0.89, 2.40)
  bResidenceTime ~ dunif(7, 14)

  sAbundanceYear ~ dunif(0, 5)
  sPeakTiming ~ dunif(0, 14)
  for(i in 1:nYear){
    bAbundanceYear[i] ~ dnorm(0, sAbundanceYear^-2)
    bPeakTimingYear[i] ~ dnorm(0, sPeakTiming^-2)
  }

  bSD ~ dnorm(0, 10^-2)
  bSDAbundance ~ dnorm(0, 5^-2)
  for(i in 1:length(Count)){

    log(eAbundance[i]) <- bAbundance + bAbundanceYear[Year[i]]
    ePeakTiming[i] <- bPeakTiming + bPeakTimingYear[Year[i]]
    log(eDuration[i]) <- bDuration

    eArrived[i] <- pnorm(Dayte[i], (ePeakTiming[i] - bResidenceTime/2), eDuration[i]^-2)
    eDeparted[i] <- pnorm(Dayte[i], (ePeakTiming[i] + bResidenceTime/2), eDuration[i]^-2)
    eSpawners[i] <- (eArrived[i] - eDeparted[i]) * eAbundance[i]

    eEfficiency[i] <- bEfficiency

    log(eSD[i]) <- bSD + bSDAbundance * eAbundance[i]

    Count[i] ~ dnorm(eSpawners[i] * eEfficiency[i], eSD[i]^-2)
  }
..

Template 2. Model definition.