Mackenzie Wood Bison Population Dynamics Analysis 2015

Suggested Citation: Thorley, J.L. and Boulanger, J. (2015) Mackenzie Wood Bison Population Dynamics Analysis 2015. A Poisson Consulting Analysis Report. URL: http://www.poissonconsulting.ca/f/164478860.

Background

The Mackenzie Wood Bison (Bison bison athabascae) herd abundance has been estimated in four years since 1999 while herd composition data has been collected in all but three years. The herd composition data is collected in July while the abundance estimates are for March.

The primary questions this analysis attempts to answer are:

What is the survival of calves, yearlings and adult in the Mackenzie herd?

Is survival of calves in the Mackenzie herd driven by climatic conditions?

Methods

Data Preparation

The data were provided by Aurora Consulting, the University of Alberta and the Government of the Northwest Territories.

In 2013 the herd experienced high mortality due to an anthrax outbreak. Consequently, the 2013 and 2014 data were excluded.

Statistical Analysis

Hierarchical Bayesian models were fitted to the data using R version 3.2.0 (Team 2013) and JAGS 3.4.0 (Plummer 2012) which interfaced with each other via jaggernaut 2.2.10 (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 specified, the models assumed vague (low information) prior 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). 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 insignificant (Kery and Schaub 2011, 37, 42) fixed (Kery and Schaub 2011, 77–82) variables and uninformative random variables. A fixed variables was considered to be insignificant if its significance was \(\geq\) 0.05 while a random variable was considered to be uninformative if its percent relative error was \(\geq\) 80%.

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 Description

The following model description observes several conventions.

The survival between life stages was estimated using a hierarchical Bayesian population dynamic state-space model (Kery and Schaub 2011). Key assumptions of the population dynamic model include:

In addition the effect of various environmental variables on calf survival was tested by adding a standardised covariate to the population dynamic model.

Model Code

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

Population Dynamic

Variable/Parameter Description
bAdults[i] Number of adults at the start of the i\(^{th}\) year
bAdults1 Number of adults at the start of the first year
bCalves[i] Number of calves at the start of the i\(^{th}\) year
Bison[i] The i\(^{th}\) herd size estimate
bProductivity Probability of a female adult calving
bSurvivalAdult Adult and yearling survival
bSurvivalCalf Calf survival
bYearlings[i] Number of yearlings at the start of the i\(^{th}\) year
bYearlings1 Number of yearlings at the start of the first year
Calves[i] Number of calves in the i\(^{th}\) composition observation
Cows[i] Number of cows in the i\(^{th}\) composition observation
Dayte[i] Day of the year of the i\(^{th}\) composition observation
eCorrection Survival correction for the timing of the herd size estimates
eProportionCalves[i] Expected proportion of cows with a calf in the i\(^{th}\) composition observation
eProportionCowsYearlings[i] Expected proportion of cows and yearlings that are cows in the i\(^{th}\) composition observation
eSurvivalCalfYear[i] Calf survival from the i\(^{th}\) to i+1\(^{th}\) year
sDispersionCalves SD of the extra-binomial variation in cow with calf clustering
sSurvivalCalfYear SD of the effect of year on bSurvivalCalf
YearBison[i] The year of the i\(^{th}\) herd size estimate
YearlingsCows[i] Number of yearlings and cows in the i\(^{th}\) composition observation
Population Dynamic - Model1
model{

  bProductivity ~ dunif(0, 1)
  bSurvivalAdult ~ dunif(0, 1)
  bSurvivalCalf ~ dunif(0, 1)

  sSurvivalCalfYear ~ dunif(0, 2)
  for(i in 1:nYear){
    bSurvivalCalfYear[i] ~ dnorm(0, sSurvivalCalfYear^-2)
    logit(eSurvivalCalfYear[i]) <- logit(bSurvivalCalf) + bSurvivalCalfYear[i]
  }

  bYearlings1 ~ dunif(0, 500)
  bAdults1 ~ dunif(0, 4000)

  bCalves[1] <- bAdults1 / 2 * bProductivity
  bYearlings[1] <- bYearlings1
  bAdults[1] <- bAdults1

  for(i in 2:nYear){
    bCalves[i] <- bAdults[i-1] / 2 * bSurvivalAdult * bProductivity
    bYearlings[i] <- bCalves[i-1] * eSurvivalCalfYear[i-1]
    bAdults[i] <- (bYearlings[i-1] + bAdults[i-1]) * bSurvivalAdult
  }

  eCorrection <- 308/365
  for(i in 1:length(YearBison)) {
    eCalves[i] <- bCalves[YearBison[i]] * eSurvivalCalfYear[YearBison[i]]^eCorrection
    eYearlings[i] <- bYearlings[YearBison[i]] * bSurvivalAdult^eCorrection
    eAdults[i] <- bAdults[YearBison[i]] * bSurvivalAdult^eCorrection

    eBison[i] <- eCalves[i] + eYearlings[i] + eAdults[i]
    Bison[i] ~ dnorm(eBison[i], 250^-2)
  }

  sDispersionCalves ~ dunif(0, 2)
  for(i in 1:length(Year)) {
    eCorComp[i] <- ((Dayte[i] - 135) / 365)
    eCalvesComp[i] <- bCalves[Year[i]] * eSurvivalCalfYear[Year[i]]^eCorComp[i]
    eYearlingsComp[i] <- bYearlings[Year[i]] * bSurvivalAdult^eCorComp[i]
    eAdultsComp[i] <- bAdults[Year[i]] * bSurvivalAdult^eCorComp[i]

    eCowsComp[i] <- eAdultsComp[i] / 2

    eDispersionCalves[i] ~ dnorm(0, sDispersionCalves^-2)
    logit(eProportionCalves[i]) <- logit(eCalvesComp[i] / eCowsComp[i]) + eDispersionCalves[i]
    eProportionCowsYearlings[i] <- eCowsComp[i] / (eYearlingsComp[i] + eCowsComp[i])

    Calves[i] ~ dbin(eProportionCalves[i], Cows[i])
    Cows[i] ~ dbin(eProportionCowsYearlings[i], YearlingsCows[i])
  }
}

Environmental

Variable/Parameter Description
bSurvivalCalfEnv Effect of the environmental variable on bSurvivalCalf
Environmental - Model1
model{

  bProductivity ~ dunif(0, 1)
  bSurvivalAdult ~ dunif(0, 1)
  bSurvivalCalf ~ dunif(0, 1)

  bSurvivalCalfEnv ~ dnorm(0, 2^-2)
  sSurvivalCalfYear ~ dunif(0, 2)
  for(i in 1:nYear){
    bSurvivalCalfYear[i] ~ dnorm(0, sSurvivalCalfYear^-2)
    logit(eSurvivalCalfYear[i]) <- logit(bSurvivalCalf) + bSurvivalCalfEnv * Env[i] + bSurvivalCalfYear[i]
  }

  bYearlings1 ~ dunif(0, 500)
  bAdults1 ~ dunif(0, 4000)

  bCalves[1] <- bAdults1 / 2 * bProductivity
  bYearlings[1] <- bYearlings1
  bAdults[1] <- bAdults1

  for(i in 2:nYear){
    bCalves[i] <- bAdults[i-1] / 2 * bSurvivalAdult * bProductivity
    bYearlings[i] <- bCalves[i-1] * eSurvivalCalfYear[i-1]
    bAdults[i] <- (bYearlings[i-1] + bAdults[i-1]) * bSurvivalAdult
  }

  eCorrection <- 308/365
  for(i in 1:length(YearBison)) {
    eCalves[i] <- bCalves[YearBison[i]] * eSurvivalCalfYear[YearBison[i]]^eCorrection
    eYearlings[i] <- bYearlings[YearBison[i]] * bSurvivalAdult^eCorrection
    eAdults[i] <- bAdults[YearBison[i]] * bSurvivalAdult^eCorrection

    eBison[i] <- eCalves[i] + eYearlings[i] + eAdults[i]
    Bison[i] ~ dnorm(eBison[i], 250^-2)
  }

  sDispersionCalves ~ dunif(0, 2)
  for(i in 1:length(Year)) {
    eCorComp[i] <- ((Dayte[i] - 135) / 365)
    eCalvesComp[i] <- bCalves[Year[i]] * eSurvivalCalfYear[Year[i]]^eCorComp[i]
    eYearlingsComp[i] <- bYearlings[Year[i]] * bSurvivalAdult^eCorComp[i]
    eAdultsComp[i] <- bAdults[Year[i]] * bSurvivalAdult^eCorComp[i]

    eCowsComp[i] <- eAdultsComp[i] / 2

    eDispersionCalves[i] ~ dnorm(0, sDispersionCalves^-2)
    logit(eProportionCalves[i]) <- logit(eCalvesComp[i] / eCowsComp[i]) + eDispersionCalves[i]
    eProportionCowsYearlings[i] <- eCowsComp[i] / (eYearlingsComp[i] + eCowsComp[i])

    Calves[i] ~ dbin(eProportionCalves[i], Cows[i])
    Cows[i] ~ dbin(eProportionCowsYearlings[i], YearlingsCows[i])
  }
}

Results

Model Parameters

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

Population Dynamic

Parameter Estimate Lower Upper SD Error Significance
bAdults1 1780.00000 1273.00000 2303.00000 259.00000 29 0.001
bProductivity 0.45970 0.40530 0.51670 0.02960 12 0.001
bSurvivalAdult 0.88991 0.85823 0.92207 0.01622 4 0.001
bSurvivalCalf 0.44130 0.32300 0.60060 0.07020 31 0.001
bYearlings1 279.70000 185.10000 381.60000 48.40000 35 0.001
sDispersionCalves 0.74820 0.60320 0.91880 0.08010 21 0.001
sSurvivalCalfYear 0.66580 0.33460 1.23600 0.23260 68 0.001
Convergence Iterations
1.01 1e+05

Environmental - Pacific Decadal Oscillation

Parameter Estimate Lower Upper SD Error Significance
bAdults1 1763.0000 1283.00000 2290.00000 261.00000 29 0.0010
bProductivity 0.4550 0.40150 0.51470 0.02940 12 0.0010
bSurvivalAdult 0.8910 0.86025 0.92222 0.01589 3 0.0010
bSurvivalCalf 0.4459 0.33070 0.57910 0.06430 28 0.0010
bSurvivalCalfEnv -0.3077 -0.73910 0.13260 0.22140 140 0.1378
bYearlings1 277.6000 191.50000 379.90000 47.70000 34 0.0010
sDispersionCalves 0.7528 0.60290 0.93010 0.08200 22 0.0010
sSurvivalCalfYear 0.6090 0.24400 1.28700 0.25900 86 0.0010
Convergence Iterations
1.01 1e+05

Environmental - Winter Severity Index

Parameter Estimate Lower Upper SD Error Significance
bAdults1 1768.00000 1293.00000 2306.00000 260.00000 29 0.0010
bProductivity 0.45660 0.40040 0.51560 0.02950 13 0.0010
bSurvivalAdult 0.88925 0.85634 0.92127 0.01612 4 0.0010
bSurvivalCalf 0.44660 0.31550 0.59530 0.07000 31 0.0010
bSurvivalCalfEnv 0.17980 -0.28890 0.67070 0.24050 270 0.3893
bYearlings1 276.70000 189.90000 374.20000 48.80000 33 0.0010
sDispersionCalves 0.74890 0.59600 0.91360 0.08140 21 0.0010
sSurvivalCalfYear 0.69000 0.32500 1.34500 0.25000 74 0.0010
Convergence Iterations
1.01 1e+05

Environmental - Rainfall

Parameter Estimate Lower Upper SD Error Significance
bAdults1 1785.50000 1294.10000 2285.60000 256.10000 28 0.0010
bProductivity 0.46020 0.40370 0.52020 0.02970 13 0.0010
bSurvivalAdult 0.88934 0.85726 0.92169 0.01665 4 0.0010
bSurvivalCalf 0.44340 0.30310 0.61400 0.07660 35 0.0010
bSurvivalCalfEnv -0.00900 -0.60000 0.56900 0.29700 6700 0.9781
bYearlings1 281.70000 189.90000 387.00000 50.00000 35 0.0010
sDispersionCalves 0.74970 0.59400 0.92190 0.07990 22 0.0010
sSurvivalCalfYear 0.71800 0.35500 1.37300 0.26000 71 0.0010
Convergence Iterations
1.01 1e+05

Environmental - Summer Air Temperature

Parameter Estimate Lower Upper SD Error Significance
bAdults1 1771.20000 1301.40000 2269.50000 252.70000 27 0.0010
bProductivity 0.45900 0.40270 0.51810 0.02970 13 0.0010
bSurvivalAdult 0.89212 0.86129 0.92205 0.01558 3 0.0010
bSurvivalCalf 0.42610 0.30690 0.57820 0.06710 32 0.0010
bSurvivalCalfEnv 0.31410 -0.12770 0.75160 0.21490 140 0.1378
bYearlings1 277.30000 197.10000 377.40000 46.40000 33 0.0010
sDispersionCalves 0.74130 0.59450 0.91120 0.08130 21 0.0010
sSurvivalCalfYear 0.60890 0.28370 1.12560 0.22700 69 0.0010
Convergence Iterations
1.01 1e+05

Figures

Population Dynamic

figures/bison/calf_cow.png
Figure 1. The predicted calf:cow ratio and composition data by year.
figures/bison/yearling_cow.png
Figure 2. The predicted yearling:cow ratio by year.
figures/bison/herd.png
Figure 3. The predicted herd size in March by year.
figures/bison/scalf.png
Figure 4. The predicted calf survival by year.

Environmental

figures/env/senv.png
Figure 5. Standardised environmental variables by year.

Acknowledgements

The organisations and individuals whose contributions have made this analysis report possible include:

References

Bradford, Michael J, Josh Korman, and Paul S Higgins. 2005. “Using Confidence Intervals to Estimate the Response of Salmon Populations (Oncorhynchus Spp.) to Experimental Habitat Alterations.” Canadian Journal of Fisheries and Aquatic Sciences 62 (12): 2716–26. https://doi.org/10.1139/f05-179.

Kery, Marc, and Michael Schaub. 2011. Bayesian Population Analysis Using WinBUGS : A Hierarchical Perspective. Boston: Academic Press. http://www.vogelwarte.ch/bpa.html.

Plummer, Martyn. 2012. “JAGS Version 3.3.0 User Manual.” http://sourceforge.net/projects/mcmc-jags/files/Manuals/3.x/.

Team, R Core. 2013. “R: A Language and Environment for Statistical Computing.” Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org.

Thorley, J. L. 2013. “Jaggernaut: An R Package to Facilitate Bayesian Analyses Using JAGS (Just Another Gibbs Sampler).” Nelson, B.C.: Poisson Consulting Ltd. https://github.com/poissonconsulting/jaggernaut.