Dewatering

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