Middle Columbia River Fish Indexing Analysis 2014

Life-Stage

The four primary fish species were categorized as fry, juvenile or adult based on their lengths.

Growth

Annual growth was estimated from the inter-annual recaptures using the Fabens method (Fabens 1965) for estimating the von Bertalanffy (VB) growth curve (von Bertalanffy 1938). The VB curves is based on the premise that

\[ \frac{dl}{dt} = k (L_{\infty} - l)\]

where \(l\) is the length of the individual, \(k\) is the growth coefficient and \(L_{\infty}\) is the mean maximum length.

Integrating the above equation gives

\[ l_t = L_{\infty} (1 - e^{-k(t - t0)})\]

where \(l_t\) is the length at time \(t\) and \(t0\) is the time at which the individual would have had no length.

The Fabens form allows

\[ l_r = L_c + (L_{\infty} - L_c) (1 - e^{-kT})\]

where \(l_r\) is the length at recapture, \(l_c\) is the length at capture and \(T\) is the time at large.

Key assumptions of the growth model include:

• \(L_{\infty}\) is constant.
• \(k\) can vary with discharge regime.
• \(k\) can vary randomly with year.
• The residual variation in growth is independently and identically normally distributed.

Condition

Condition was estimated via an analysis of mass-length relations (He et al. 2008).

More specifically the model was based on the allometric relationship

\[ W = \alpha L^{\beta}\]

where \(W\) is the weight (mass), \(\alpha\) is the coefficent, \(\beta\) is the exponent and \(L\) is the length.

To improve chain mixing the relation was log-transformed, i.e.,

\[ log(W) = log(\alpha) + \beta * log(L)\]

and the logged lengths centered, i.e., \(log(L) - \bar{log(L)}\), prior to model fitting.

Preliminary analyses indicated that the variation in the exponent \(\beta\) with respect to year was not informative.

Key assumptions of the final condition model include:

• The expected weight varies with length as an allometric relationship.
• The intercept of the log-transformed allometric relationship is described by a linear mixed model.
• The intercept of the log-transformed allometric relationship varies with discharge regime and season.
• The intercept of the log-transformed allometric relationship varies randomly with year, site and the interaction between year and site.
• The slope of the log-transformed allometric relationship is described by a linear mixed model.
• The slope of the log-transformed allometric relationship varies with discharge regime and season.
• The slope of the log-transformed allometric relationship varies randomly with year.
• The residual variation in weight for the log-transformed allometric relationship is independently and identically normally distributed.

Occupancy

Occupancy, which is the probability that a particular species was present at a site, was estimated from the temporal replication of detection data (Kery, 2010; Kery and Schaub, 2011, pp. 238-242 and 414-418), i.e., each site was surveyed multiple times within a season. A species was considered to have been detected if one or more individuals of the species were caught or counted. It is important to note that the model estimates the probability that the species was present at a given (or typical) site in a given (or typical) year as opposed to the probability that the species was present in the entire study area. We focused on Northern Pikeminnow, Burbot, Lake Whitefish, Rainbow Trout, Redside Shiner and Sculpins because they were low enough density to not to be present at all sites at all times yet were encounted sufficiently often to provide information on spatial and temporal changes.

Key assumptions of the occupancy model include:

• Occupancy (probability of presence) is described by a generalized linear mixed model with a logit link.
• Occupancy varies with discharge regime and season.
• Occupancy varies randomly with site and year.
• Sites are closed, i.e., the species is present or absent at a site for all the sessions in a particular season of a year.
• Observed presence is described by a bernoulli distribution, given occupancy.

Species Richness

The estimated probabilities of presence for the six species considered in the occupany analyses were summed to give the expected species richnesses at a given (or typical) site in a given (or typical) year.

Count

The count data were analysed using an overdispersed Poisson model (Kery, 2010; Kery and Schaub, 2011, pp. 168-170,180 and 55-56) to provide estimates of the lineal river count density (count/km) by year and site. Unlike Kery (2010) and Kery and Schaub (2011), which used a log-normal distribution to account for the extra-Poisson variation, the current model used a gamma distribution with identical shape and scale parameters because it has a mean of 1 and therefore no overall effect on the expected count. The count data does not enable us to estimate abundance nor observer efficiency, but it enables us to estimate an expected count, which is the product of the two. As such it is necessary to assume that changes in observer efficiency are negligible in order to interpret the estimates as relative density.

Key assumptions of the abundance model include:

• Lineal density (fish/km) is described by a generalized linear mixed model with a logarithm link.
• Lineal density (fish/km) varies with discharge regime and season.
• Lineal density (fish/km) varies randomly with river kilometer, site, year and the interaction between site and year.
• The regression coefficient of river kilometre is described by a linear mixed model.
• The effect of river kilometre on lineal density (distribution of density along the river) varies with discharge regime and season.
• The effect of river kilometre on lineal density (distribution of density along the river) varies randomly with year.
• The counts are gamma-Poisson distributed, given the mean count.

Movement

The extent to which sites are closed, i.e., fish remain at the same site between sessions, was evaluated from a logistic ANCOVA (Kery 2010). The model estimated the probability that intra-annual recaptures were caught at the same site versus a different one.

Key assumptions of the site fidelity model include:

• Site fidelity varies with season, length and the interaction between season and length.
• Observed site fidelity is Bernoulli distributed.

Observer Length Correction

The bias (accuracy) and error (precisions) in observer’s fish length estimates were quantified using a model with a categorical distribution that compared the proportions of fish in different length-classes for each observer to the equivalent proportions for the measured fish.

Key assumptions of the observer length correction model include:

• The expected length bias can vary by observer.
• The expected length error can vary by observer.
• The residual variation in length is independently and identically normally distributed.

The observed fish lengths were corrected for the estimated length biases.

Abundance

The catch data were analysed using a capture-recapture-based overdispersed Poisson model to provide estimates of capture efficiency and absolute abundance.
To maximize the number of recaptures the model grouped all the sites into a supersite for the purposes of estimating the number of marked fish but analysed the total captures at the site level.

Key assumptions of the abundance model include:

• Lineal density (fish/km) varies with discharge regime, season and river km.
• Lineal density (fish/km) varies randomly with site, year and the interaction between site and year.
• The relationship between density and river kilometre (distribution) varies with discharge regime and season.
• The relationship between density and river kilometre (distribution) varies randomly with year.
• Efficiency (probability of capture) varies by season and method (capture versus count).
• Efficiency varies randomly by session within season within year.
• Marked and unmarked fish have the same probability of capture.
• There is no tag loss, migration (at the supersite level), mortality or misidentification of fish.
• The number of fish caught is gamma-Poisson distributed.

Species Evenness

The site and year estimates of the lineal bank count densities from the count model for Rainbow Trout, Suckers, Burbot and Northern Pikeminnow were combined with the equivalent count estimates for Bull Trout and Adult Mountain Whitefish from the abundance model to calculate the shannon index of evenness (\(E\)). The index was calculated using the following formula where \(S\) is the number of species and \(p_i\) is the proportion of the total count belonging to the ith species.

\[ E = \frac{-\sum p_i \log(p_i)}{ln(S)}\]

Long-Term Trends

The long-term congruence between the yearly fall fish metrics (Grw = Growth, Con = Condition, Blank = Abundance, Rkm = Distribution) by life stage (Juv = Juvenile, AD = Adult) and a range of environmental variables including the tri-monthly (JaMa = January - March, ApJn = April - June, JlSe = July - September, OcDe = October - December) discharge from Revelstoke Dam (Q10, QMu = Mean, Q90, QDlt = Mean absolute difference) and the elevation (Ele), the Pacific Decadal Oscillation (PDO), the chlorophyll-a (ChlA) and invertebrate production (Inv) and the kokanee abundance (KO) was examined using Dynamic Factor Analysis (Zuur, Tuck, and Bailey 2003). Due to the rotation problem the underlying trends were indeterminate (Abmann, Boysen-Hogrefe, and Pape 2014). The October to December discharge and elevation values were lead (opposite of lagged) by one year to account for the fact that they occurred after sampling.

Key assumptions of the Dynamic Factor Analysis (DFA) model include:

• The time series are described by three underlying trends.
• The random walk processes in the trends are normally distributed.
• The residual variation in the standardised variables is normally distributed.

The similarities were represented visually by using non-metric multidimensional scaling (NMDS) to map the distances onto two-dimensional space. The more similar two time series are they closer they will tend to be in the resultant NMDS plot.

Short-Term Correlations

To assess the short-term congruence between the yearly fish metrics and the environmental variables, the pair-wise distances between the residuals from the DFA model were calculated as \(1 - abs(cor(x, y))\) where \(cor\) is the Pearson correlation, \(abs\) the absolute value and \(x\) and \(y\) are the two time series being compared.

The short-term similarities were represented visually by using NMDS to map the distances onto two-dimensional space.

Growth

model {

  bKIntercept ~ dnorm (0, 5^-2)

  bKRegime[1] <- 0
  for(i in 2:nRegime) {
    bKRegime[i] ~ dnorm(0, 5^-2)
  }

  sKYear ~ dunif (0, 5)
  for (i in 1:nYear) {
    bKYear[i] ~ dnorm(0, sKYear^-2)
    log(bK[i]) <- bKIntercept + bKRegime[step(i - Threshold) + 1] + bKYear[i]
  }

  bLinf ~ dunif(100, 1000)
  sGrowth ~ dunif(0, 100)

  for (i in 1:length(Year)) {

    eGrowth[i] <- (bLinf - LengthAtRelease[i]) * (1 - exp(-sum(bK[Year[i]:(Year[i] + Years[i] - 1)])))

    Growth[i] ~ dnorm(eGrowth[i], sGrowth^-2)
  }
  tGrowth <- bKRegime[2]
} 

Condition

model {

  bWeightIntercept ~ dnorm(5, 5^-2)
  bWeightSlope ~ dnorm(3, 5^-2)

  bWeightRegimeIntercept[1] <- 0
  bWeightRegimeSlope[1] <- 0

  for(i in 2:nRegime) {
    bWeightRegimeIntercept[i] ~ dnorm(0, 5^-2)
    bWeightRegimeSlope[i] ~ dnorm(0, 5^-2)
  }

  bWeightSeasonIntercept[1] <- 0
  bWeightSeasonSlope[1] <- 0
  for(i in 2:nSeason) {
    bWeightSeasonIntercept[i] ~ dnorm(0, 5^-2)
    bWeightSeasonSlope[i] ~ dnorm(0, 5^-2)
  }

  sWeightYearIntercept ~ dunif(0, 5)
  sWeightYearSlope ~ dunif(0, 5)
  for(yr in 1:nYear) {
    bWeightYearIntercept[yr] ~ dnorm(0, sWeightYearIntercept^-2)
    bWeightYearSlope[yr] ~ dnorm(0, sWeightYearSlope^-2)
  }

  sWeightSiteIntercept ~ dunif(0, 5)
  sWeightSiteYearIntercept ~ dunif(0, 5)
  for(st in 1:nSite) {
    bWeightSiteIntercept[st] ~ dnorm(0, sWeightSiteIntercept^-2)
    for(yr in 1:nYear) {
      bWeightSiteYearIntercept[st, yr] ~ dnorm(0, sWeightSiteYearIntercept^-2)
    }
  }

  sWeight ~ dunif(0, 5)
  for(i in 1:length(Year)) {

    eWeightIntercept[i] <- bWeightIntercept
        + bWeightRegimeIntercept[Regime[i]]
        + bWeightSeasonIntercept[Season[i]]
        + bWeightYearIntercept[Year[i]]
        + bWeightSiteIntercept[Site[i]]
        + bWeightSiteYearIntercept[Site[i],Year[i]]

    eWeightSlope[i] <- bWeightSlope
        + bWeightRegimeSlope[Regime[i]]
        + bWeightSeasonSlope[Season[i]]
        + bWeightYearSlope[Year[i]]

    log(eWeight[i]) <- eWeightIntercept[i] + eWeightSlope[i] * Length[i]
    Weight[i] ~ dlnorm(log(eWeight[i]) , sWeight^-2)
  }
  tCondition1 <- bWeightRegimeIntercept[2]
  tCondition2 <- bWeightRegimeSlope[2]
}

Occupancy

model {

  bOccupancy ~ dnorm(0, 5^-2)
  bOccupancySeason[1] <- 0
  for(i in 2:nSeason) {
    bOccupancySeason[i] ~ dnorm(0, 5^-2)
  }

  bOccupancyRegime[1] <- 0
  for(i in 2:nRegime) {
    bOccupancyRegime[i] ~ dnorm(0, 5^-2)
  }

  sOccupancyYear ~ dunif(0, 5)
  for (yr in 1:nYear) {
    bOccupancyYear[yr] ~ dnorm(0, sOccupancyYear^-2)
  }

  sOccupancySite ~ dunif(0, 5)
  for (st in 1:nSite) {
    bOccupancySite[st] ~ dnorm(0, sOccupancySite^-2)
  }
  for (i in 1:length(Year)) {
    logit(eOccupancy[i]) <- bOccupancy
    + bOccupancyRegime[Regime[i]] + bOccupancySeason[Season[i]]
    + bOccupancySite[Site[i]] + bOccupancyYear[Year[i]]
    eObserved[i] <- eOccupancy[i]
    Observed[i] ~ dbern(eObserved[i])
  }
}

Count

model {
  bDensity ~ dnorm(0, 5^-2)
  bDistribution ~ dnorm(0, 5^-2)
  bDensityRegime[1] <- 0

  bDistributionRegime[1] <- 0
  for(i in 2:nRegime) {
    bDensityRegime[i] ~ dnorm(0, 5^-2)
    bDistributionRegime[i] ~ dnorm(0, 5^-2)
  }

  bDensitySeason[1] <- 0
  bDistributionSeason[1] <- 0
  for(i in 2:nSeason) {
    bDensitySeason[i] ~ dnorm(0, 5^-2)
    bDistributionSeason[i] ~ dnorm(0, 5^-2)
  }

  sDensityYear ~ dunif(0, 2)
  sDistributionYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bDensityYear[i] ~ dnorm(0, sDensityYear^-2)
    bDistributionYear[i] ~ dnorm(0, sDistributionYear^-2)
  }

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

  sDispersion ~ dunif(0, 5)
  for (i in 1:length(Year)) {
    eDistribution[i] <- bDistribution
      + bDistributionRegime[Regime[i]]
      + bDistributionSeason[Season[i]]
      + bDistributionYear[Year[i]]

    log(eDensity[i]) <- bDensity
      + eDistribution[i] * RiverKm[i]
      + bDensityRegime[Regime[i]]
      + bDensitySeason[Season[i]]
      + bDensitySite[Site[i]]
      + bDensityYear[Year[i]]
      + bDensitySiteYear[Site[i],Year[i]]

    eCount[i] <- eDensity[i] * SiteLength[i] * ProportionSampled[i]
    eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)
    Count[i] ~ dpois(eCount[i] * eDispersion[i])
  }
  tAbundance <- bDensityRegime[2]
  tDistribution <- bDistributionRegime[2]
}

Movement

model {
  bMoved ~ dnorm(0, 5^-2)
  bLength ~ dnorm(0, 5^-2)
  bMovedSeason[1] <- 0
  bLengthSeason[1] <- 0

  for(i in 2:nSeason) {
    bMovedSeason[i] ~ dnorm(0, 5^-5)
    bLengthSeason[i] ~ dnorm(0, 5^-5)
  }

  for (i in 1:length(Season)) {
    logit(eMoved[i]) <- bMoved + bMovedSeason[Season[i]] + (bLength + bLengthSeason[Season[i]]) * Length[i]
    Moved[i] ~ dbern(eMoved[i])
  }
}

Observer Length Correction

model {
  for(i in 1:nClass) {
    dClass[i] <- 1
  }
  pClass[1:nClass] ~ ddirch(dClass[])

  bLength[1] <- 1
  sLength[1] <- 1

  for(i in 2:nObserver) {
    bLength[i] ~ dunif(0.5, 2)
    sLength[i] ~ dunif(1, 50)
  }
  for(i in 1:length(Length))  {
    eClass[i] ~ dcat(pClass[])
    eLength[i] <- bLength[Observer[i]] * ClassLength[eClass[i]]
    eSLength[i] <- sLength[Observer[i]] * ClassSD
    Length[i] ~ dnorm(eLength[i], eSLength[i]^-2)
  }
}

Abundance

model {

  bEfficiency ~ dnorm(0, 5^-2)
  bDensity ~ dnorm(0, 5^-2)
  bDistribution ~ dnorm(0, 5^-2)

  bDensityRegime[1] <- 0
  bDistributionRegime[1] <- 0
  for(i in 2:nRegime) {
    bDensityRegime[i] ~ dnorm(0, 5^-2)
    bDistributionRegime[i] ~ dnorm(0, 5^-2)
  }

  bEfficiencySeason[1] <- 0
  bDensitySeason[1] <- 0
  bDistributionSeason[1] <- 0
  for(i in 2:nSeason) {
    bEfficiencySeason[i] ~ dnorm(0, 5^-2)
    bDensitySeason[i] ~ dnorm(0, 5^-2)
    bDistributionSeason[i] ~ dnorm(0, 5^-2)
  }

  sDensityYear ~ dunif(0, 2)
  sDistributionYear ~ dunif(0, 2)
  for (i in 1:nYear) {
    bDensityYear[i] ~ dnorm(0, sDensityYear^-2)
    bDistributionYear[i] ~ dnorm(0, sDistributionYear^-2)
  }

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

  sEfficiencySessionSeasonYear ~ dunif(0, 5)
  for (i in 1:nSession) {
    for (j in 1:nSeason) {
      for (k in 1:nYear) {
        bEfficiencySessionSeasonYear[i, j, k] ~ dnorm(0, sEfficiencySessionSeasonYear^-2)
      }
    }
  }

  bType[1] <- 1
  for (i in 2:nType) {
    bType[i] ~ dunif(0, 10)
  }

  for(i in 1:length(EffIndex)) {

    logit(eEff[i]) <- bEfficiency
        + bEfficiencySeason[Season[EffIndex[i]]]
        + bEfficiencySessionSeasonYear[Session[EffIndex[i]],
                                       Season[EffIndex[i]],
                                       Year[EffIndex[i]]]

    Marked[EffIndex[i]] ~ dbin(eEff[i], Tagged[EffIndex[i]])
  }

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

    logit(eEfficiency[i]) <- bEfficiency
        + bEfficiencySeason[Season[i]]
        + bEfficiencySessionSeasonYear[Session[i], Season[i], Year[i]]

    eDistribution[i] <- bDistribution
      + bDistributionRegime[Regime[i]]
      + bDistributionSeason[Season[i]]
      + bDistributionYear[Year[i]]

    log(eDensity[i]) <- bDensity
      + eDistribution[i] * RiverKm[i]
      + bDensityRegime[Regime[i]]
      + bDensitySeason[Season[i]]
      + bDensitySite[Site[i]]
      + bDensityYear[Year[i]]
      + bDensitySiteYear[Site[i], Year[i]]

    eCatch[i] <- eDensity[i] * SiteLength[i] * ProportionSampled[i] * eEfficiency[i] * bType[Type[i]]

    eDispersion[i] ~ dgamma(1 / sDispersion^2, 1 / sDispersion^2)

    Catch[i] ~ dpois(eCatch[i] * eDispersion[i])
  }
  tAbundance <- bDensityRegime[2]
  tDistribution <- bDistributionRegime[2]
}

Long-Term Trends

model{

  sTrend ~ dunif(0, 1)
  for (t in 1:nTrend) {
    bTrendYear[t,1] ~ dunif(-1,1)
    for(y in 2:nYear){
      bTrendYear[t,y] ~ dnorm(bTrendYear[t,y-1], sTrend^-2)
    }
  }

  for(v in 1:nVariable){
    for(t in 1:nTrend) {
      Z[v,t] ~ dunif(-1,1)
    }
    for(y in 1:nYear){
      eValue[v,y,1] <- Z[v,1] * bTrendYear[1,y]
      for(t in 2:nTrend) {
        eValue[v,y,t] <- eValue[v,y,t-1] + Z[v,t] * bTrendYear[t,y]
      }
    }
  }

  sValue ~ dunif(0, 1)
  for(i in 1:length(Value)) {
    Value[i] ~ dnorm(eValue[Variable[i], Year[i], nTrend], sValue^-2)
  }

  for(i in 1:nVariable) {
    for(j in 1:nVariable) {
      bDistance[i,j] <- sqrt(sum((Z[i,]-Z[j,])^2))
    }
  }
}

Growth - Bull Trout

Growth - Mountain Whitefish

Growth - Rainbow Trout

Condition - Bull Trout

Condition - Mountain Whitefish

Condition - Rainbow Trout

Occupancy - Rainbow Trout

Occupancy - Burbot

Occupancy - Lake Whitefish

Occupancy - Northern Pikeminnow

Occupancy - Redside Shiner

Occupancy - Sculpins

Count - Rainbow Trout

Count - Burbot

Count - Northern Pikeminnow

Count - Suckers

Movement - Bull Trout

Movement - Mountain Whitefish

Movement - Rainbow Trout

Movement - Largescale Sucker

Observer Length Correction - Bull Trout

Observer Length Correction - Mountain Whitefish

Observer Length Correction - Suckers

Abundance - Bull Trout - Adult

Abundance - Bull Trout - Juvenile

Abundance - Mountain Whitefish - Adult

Abundance - Mountain Whitefish - Juvenile

Abundance - Rainbow Trout - Adult

Abundance - Largescale Sucker - Adult

Long-Term Trends