Lower Columbia River Rainbow Trout Spawning Analysis 2016

Data Preparation

The Rainbow Trout fish and redd aerial count data for the LCR and LKR were collected by Mountain Water Research and databased by Gary Pavan. The age-1 Rainbow Trout abundance data were provided by the Fish Population Indexing Program which is led by Okanagan Nation Alliance.

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). A decline in the redd count of more than one third of the previous maximum count for a particular section was inferred to be caused by poor viewing conditions (turbidity) and the affected spawner and redd section counts were excluded from any subsequent analyses.

The data were prepared for analysis using R version 3.3.2 (R Core Team 2015).

Statistical Analysis

Hierarchical Bayesian models were fitted to the data using R version 3.3.2 (R Core Team 2015) and JAGS 4.0.1 (Plummer 2015) which interfaced with each other via jaggernaut 2.5.1 (Thorley 2013). For additional information on hierarchical Bayesian modelling in the BUGS language, of which JAGS uses a dialect, the reader is referred to Kery and Schaub (2011, 41–44).

Unless indicated otherwise, the models used prior distributions that were vague in the sense that they did not affect the posterior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from a minimum of 1,000 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 Rhat (Kery and Schaub 2011, 40) was less than 1.1 for each of the parameters in the model (Kery and Schaub 2011, 61). Where relevant, model adequacy was confirmed by examination of residual plots.

The posterior distributions of the fixed (Kery and Schaub 2011, 75) parameters are summarised in terms of a point estimate (mean), lower and upper 95% credible limits (2.5th and 97.5th percentiles), the standard deviation (SD), percent relative error (half the 95% credible interval as a percent of the point estimate) and significance (Kery and Schaub 2011, 37, 42).

Variable selection was achieved by dropping fixed (Kery and Schaub 2011, 77–82) variables with two-sided p-values \(\geq\) 0.05 (Kery and Schaub 2011, 37, 42) and random variables with percent relative errors \(\geq\) 80%. The Deviance Information Criterion (DIC) was not used because it is of questionable validity when applied to hierarchical models (Kery and Schaub 2011, 469).

The results are displayed graphically by plotting the modeled relationships between particular variables and the response with 95% credible intervals (CRIs) 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). Where informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% CRIs (Bradford, Korman, and Higgins 2005).

Model Descriptions

Stock-Recruitment

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

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

where \(S\) is the adults (stock), \(R\) is the subadults (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 prior probability distribution for the maximum number of recruits per spawner \(\alpha\) is normally distributed with a mean of 90 and a SD of 50.
• The density-dependence varies with the proportional egg loss.
• The residual variation in the number of age-1 recruits is log-normally distributed.

The prior probability distribution 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 \(K\) is given by the relationship:

\[ K = \frac{\alpha}{\beta} \quad.\]

Model Code

The JAGS model code, which uses a series of naming conventions, is presented below.

model{
  bSpAbundance ~ dnorm(5, 5^-2)

  bSpArrivalPeak ~ dnorm(0, 14^-2)
  sSpArrivalWidth ~ dunif(log(14), log(42))
  bSpResidence ~ dnorm(11, 3.07^-2) T(6.31, 18.16)

  bSpObsEfficiency ~ dunif(0.9, 1.1)

  bSpAbundanceSite[1] <- 0
  for (i in 2:nSite) {
    bSpAbundanceSite[i] ~ dnorm(0, 2^-2)
  }

  sSpAbundanceYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bSpAbundanceYear[i] ~ dnorm(0, sSpAbundanceYear^-2)
  }

  sSpAbundanceSiteYear ~ dunif(0, 2)
  for (i in 1:nSite) {
    for (j in 1:nYear) {
      bSpAbundanceSiteYear[i, j] ~ dnorm(0, sSpAbundanceSiteYear^-2)
    }
  }

  sSpArrivalPeakYear ~ dunif(0, 28)
  for (i in 1:nYear) {
    bSpArrivalPeakYear[i] ~ dnorm(0, sSpArrivalPeakYear^-2)
  }

  bSpArrivalWidthSite[1] <- 0
  for(i in 2:nSite){
    bSpArrivalWidthSite[i] ~ dnorm(0, 1^-2)
  }

  bReddPerSpawner ~ dunif(0, 4)

  bRdResidence ~ dnorm(100, 50^-2)

  bRdObsEfficiency ~ dunif(0.9, 1.1)

  sFishDispersion ~ dunif(0, 2)
  sReddDispersion ~ dunif(0, 2)

  for (i in 1:length(Fish)) {
    log(eSpAbundance[i]) <- bSpAbundance +
                            bSpAbundanceSite[Site[i]] +
                            bSpAbundanceYear[Year[i]] +
                            bSpAbundanceSiteYear[Site[i], Year[i]]

    eSpArrivalPeak[i] <- bSpArrivalPeak + bSpArrivalPeakYear[Year[i]]
    log(eSpArrivalWidth[i]) <- sSpArrivalWidth + bSpArrivalWidthSite[Site[i]]

    eSpFracArrived[i] <- pnorm(
                           Dayte[i],
                           (eSpArrivalPeak[i] - bSpResidence/2),
                           eSpArrivalWidth[i]^-2
                         )
    eSpFracDeparted[i] <- pnorm(
                            Dayte[i],
                            (eSpArrivalPeak[i] + bSpResidence/2),
                            eSpArrivalWidth[i]^-2
                          )
    eFish[i] <- (eSpFracArrived[i] - eSpFracDeparted[i])
                * eSpAbundance[i]
                * bSpObsEfficiency

    eFishDispersion[i] ~ dgamma(1/sFishDispersion^2, 1/sFishDispersion^2)
    Fish[i] ~ dpois(eFish[i] * eFishDispersion[i])

    eRdAbundance[i] <- eSpAbundance[i] * bReddPerSpawner

    eRdFracArrived[i] <- pnorm(
                           Dayte[i],
                           (eSpArrivalPeak[i] - bSpResidence/2),
                           eSpArrivalWidth[i]^-2
                         )
    eRdFracDeparted[i] <- pnorm(
                            Dayte[i],
                            (eSpArrivalPeak[i] + bRdResidence/2),
                            eSpArrivalWidth[i]^-2
                          )
    eRedds[i] <- (eRdFracArrived[i] - eRdFracDeparted[i])
                 * eRdAbundance[i]
                 * bRdObsEfficiency

    eReddDispersion[i] ~ dgamma(1/sReddDispersion^2, 1/sReddDispersion^2)
    Redds[i] ~ dpois(eRedds[i] * eReddDispersion[i])
  }
}

model{
  alpha ~ dnorm(90, 50^-2) T(1, )
  k ~  dnorm(2*10^4, (2*10^3)^-2) T(0, )
  sRecruits ~ dunif(0, 5)

  bBetaEggLoss ~ dnorm(0, 2^-2)

  for(i in 1:length(Stock)){
    log(eBeta[i]) <- log(alpha / k) + bBetaEggLoss * EggLoss[i]
    eRecruits[i] <- alpha * Stock[i] / (1 + Stock[i] * eBeta[i])
    Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
  }
}

model{
  bResidenceTime ~ dnorm(0, 5^-2)
  sResidenceTime ~ dunif(0, 5)

  for(i in 1:length(ResidenceTime)){
    log(eResidenceTime[i]) <- bResidenceTime
    ResidenceTime[i] ~ dlnorm(log(eResidenceTime[i]), sResidenceTime^-2)
  }
}

Model Parameters

The posterior distributions for the fixed (Kery and Schaub 2011 p. 75) parameters in each model are summarised below.

Figures

Discharge

Temperature

Area-Under-The-Curve

Stock-Recruitment

Acoustic