Bayesian Models for Ecologists
Dynamic Models: Forecasting The Effects of Harvest on Lynx
June 12, 2024

Motivation

The Eurasian lynx (Lynx lynx) is a medium-sized predator with broad distribution in the boreal forests of Europe and Siberia. The lynx is classified as a threatened species throughout much of its range.

There is controversy about the legal harvest of lynx in Sweden. Proponents of harvest argue that allowing hunting of lynx reduces illegal kill (poaching). Moreover, Sweden is committed to regulate lynx numbers to prevent excessive predation on reindeer because reindeer are critical to the livelihoods of indigenous pastoralists, the Sami. Many environmentalists oppose harvest, however, arguing that lynx are too rare to remove their fully protected status. A similar controversy surrounds management of wolves in the Western United States.

drawing

drawing

A forecasting model for the abundance of lynx helps managers make decisions that can be justified to citizens. The model you will develop today is not a toy. A small modfication of the model is currently used in Sweden and Norway to manage Lynx (H. Andren, N. T. Hobbs, M. Aronsson, H. Broseth, G. Chapron, J. D. C. Linnell, J. Odden, J. Persson, and E. B. Nilsen. Harvest models of small populations of a large carnivore using Bayesian forecasting. Ecological Applications, 30(3):e02063, 2020.)

You have data on the number of lynx family groups censused in a management unit as well as annual records of lynx harvested from the unit. You will model the population using the deterministic model:

\[N_t=\lambda(N_{t-1}-H_{t-1}).\]

where \(N_{t}\) is the true, unobserved abundance of lynx and \(H_{t-1}\) is the number of lynx harvested during \(t-1\) to \(t\). The parentheses in this expression reflect the fact that harvest occurs immediately after census, such that the next years population increment comes from the post-harvest population size. The term census as used here refers to the time that the population is updated in the model relative to other events in the year. The sequence of these events matters in discrete time models.

ADVANCED (for the population modelers) What would be the model if harvest occurred immediately before census? Three months after census? Continuously throughout the year?


Immediately before census:

\[N_t=\lambda N_{t-1}-H_{t-1}.\] Three months after census:

\[N_t=(\lambda^{\frac{1}{4}}N_{t-1}-H_{t-1})\lambda^{\frac{3}{4}}\]

Throughout the year:

\[N_t=(\lambda^{\frac{1}{2}}N_{t-1}-H_{t-1})\lambda^{\frac{1}{2}}\] \[N_t=\lambda N_{t-1}-\lambda^{\frac{1}{2}}H_{t-1}\]


Assume the harvest (\(H_t\)) and the number of family groups (\(y_t\) ) are unbiased such that there is no consistent under or over-counting. You are entitled to make the assumption that family groups are observed without bias because your Scandinavian colleagues are competent snow trackers and do a good job of estimating the number of family groups (if not the number of lynx) in a management region. Harvest can be assumed to be measured without error because it is closely regulated and all hunters who harvest a lynx are required by law to register the animal with the county.

The challenge in this problem is that the observations of lynx abundance (family groups) are not the same as the observation of harvest (number of lynx). Fortunately, you have prior information, hard won from radio-telemetry, on the proportional relationship between number of family groups and number of lynx in the population, i.e:

\[\phi=\frac{f}{N},\]

where \(f\) is the number of family groups and \(N\) is the population size, mean \(\phi=0.163\) with standard deviation of the mean = 0.012.


R libraries needed for this lab

You need to load the following libraries. Set the seed to 10 to compare your answers to ours. The data for this problem is located in the LynxFamilies data frame of the BayesNSF package.

library(BayesNSF)
library(rjags)
library(MCMCvis)
library(HDInterval)
set.seed(10)


Generating an Informed Prior for \(\phi\)

We’ve provided you with a useful moment matching function below for converting the mean and standard deviation of \(\phi\) to the parameters for the beta distribution you will use as an informed prior on \(\phi\).

# Function to get beta shape parameters from moments
shape_from_stats <- function(mu = mu.global, sigma = sigma.global) {
  a <-(mu^2 - mu^3 - mu * sigma^2) / sigma^2
  b <- (mu - 2 * mu^2 + mu^3 - sigma^2 + mu*sigma^2) / sigma^2
  shape_ps <- c(a, b)
  return(shape_ps)
}

# get parameters for distribution of population multiplier, 1/p
shapes = shape_from_stats(.163, .012)

# check prior on p using simulated data from beta distribution
x = seq(0, 1, .001)
p = dbeta(x, shapes[1], shapes[2])
plot(x, p, typ = "l", xlim = c(0, 1))


Diagram the Bayesian network
  1. Develop a hierarchical Bayesian model (also called a state space model) of the lynx population in the management unit. Diagram the Bayesian network (the DAG) of knowns and unknowns and write out the posterior and factored joint distribution. Use a lognormal distribution to model the true lynx population size over time. Use a Poisson distribution for the data model relating the true, unobserved state (the total population size) to the observed data (number of family groups).


\[\begin{align*} [\phi, \lambda, \pmb{N}, \sigma^{2}_{p} \mid \pmb{y}] &\propto \prod_{t=2}^{n}[y_{t} \mid N_{t}, \phi][N_{t} \mid \lambda, N_{t-1}, \sigma^{2}_{p}] [N_{1}\mid y_1, \phi, \sigma^{2}_{p}][\phi][\lambda][\sigma^{2}_{p}]\\ y_{t} &\sim \text{Poisson}(N_{t}\phi)\\ N_{t} &\sim \text{lognormal}(\log(\lambda(N_{t-1}-H_{t-1})), \sigma^{2}_{p})\\ \phi &\sim \text{beta}(154, 792)\\ N_{1} &\sim \text{lognormal}\bigg(\log\bigg(\frac{y_{1}}{\phi}\bigg), \sigma^{2}_{p}\bigg)\\ \lambda &\sim \text{uniform}(0.1, 10)\\ \sigma^{2}_{p} &\sim \text{uniform}(0, 5)\\ \end{align*}\]


  1. ADVANCED An alternative approach, which is slightly more difficult to code, is to model the process as:

\[\text{negative binomia}(N_t|\lambda(N_{t-1}-H_{t-1}, \rho)),\]

and model the data as:

\[\text{binomial}(y_t| \text{round}(N_t\phi),p),\]

where \(p\) is a detection probability. Explain why this second formulation might be better than the formulation you are using. (It turns out they give virtually identical results.)


There are two advantages to using a negative binomial process model and binomial data model. A negative binomial process model treats the true state as an integer, which for small populations, has some advantages because it includes demographic stochasticity. The binomial data model assures that the observed state is never larger than the true state, which makes sense if the only source of error in the census is failing to observe families that are present. On the other hand, if there is a possibility of double-counting, which is the case here, the then Poisson is a better choice for the data model.


Fitting the Model

Estimate the marginal posterior distribution of the unobserved, true state over time (\(\mathbf{N}\)), the parameters in the model \(\lambda\), and \(\phi\), as well as the process variance and observation variance. Summarize the marginal posterior distributions of the parameters and unobserved states.

A note about the data. Each row in the data file gives the observed number of family groups for that year in column 2 and that year’s harvest in column 3. The harvest in each row influences the population size in the next row. So, for example, the 2016 harvest influences the 2017 population size.

Before you begin it’s very helpful to use simulated data to the verify initial values and model. We simulate the true state by choosing some biologically reasonable values for model parameters and “eyeballing” the fit of the true state to the data. You can then use these simulated values for initial conditions (see the inits list below). This is of particular importance because failing to give reasonable initial conditions for dynamic models can cause problems in model fitting. Remember, supply initial conditions for all unobserved quantities in the posterior distribution (even those that do not have priors).

y <- LynxFamilies
endyr <- nrow(y)
n <- numeric(endyr + 1)
mu <- numeric(endyr + 1)
fg <- numeric(endyr + 1)
phi <- 0.16
lambda <- 1.07
sigma.p <- 0.2
 
n[1] <- y$census[1] / phi # n in the unit of individuals
mu[1] <- n[1] # mean from deterministic model to simulate
fg[1] <- n[1] * phi # Nt in the unit of
 
for (t in 2:(endyr + 1)) {
  mu[t] <- lambda * (n[t - 1] - y$harvest[t - 1])
  n[t] <- rlnorm(1, log(mu[t]), sigma.p)
  fg[t] <- n[t] * phi
}

plot(y$year, y$census, ylim = c(0, 100), xlab = "Year", ylab = "Family group", main = "Simulated data")
lines(y$year, fg[1:length(y$year)])

## visually match simulated data with observations for initial conditions
endyr = nrow(y)
n = numeric(endyr + 1)
mu = numeric(endyr + 1) #use this for family groups
lambda = 1.1
sigma.p = .00001
n[1] = y$census[1]

for(t in 2:(endyr + 1)) {
  n[t] <- lambda * (y$census[t - 1] - .16 * y$harvest[t - 1])  # use this for family groups
}

plot(y$year, y$census, ylim = c(0, 100), xlab = "Year", ylab = "Family group", main = "Simulated data")
lines(y$year, n[1:length(y$year)])

Here’s your starting code:

data = list(
    y.endyr = endyr,
    y.a = shapes[1], 
    y.b = shapes[2],
    y.H = y$harvest,
    y = y$census)

inits = list(
    list(lambda = 1.2, sigma.p = .01, N = n),
    list(lambda = 1.01,sigma.p = .2, N = n * 1.2),
    list(lambda = .95, sigma.p = .5, N = n * .5))


  1. Write the JAGS model to estimate the marginal posterior distribution of the unobserved, true state over time (\(\mathbf{N}\)), the parameters in the model \(\lambda\) and \(\phi\) as well as the process variance and observation variance. Include a summary the marginal posterior distributions of the parameters and unobserved states.


{
sink("LynxHarvest.R")
cat("

model{

  # priors
  sigma.p ~ dunif(0, 5)
  tau.p <- 1 / sigma.p^2
  lambda ~ dunif(0, 10)
  p ~ dbeta(y.a, y.b)  

  # Get parameters y.a and y.b from mean and sd using moment matching to make this prior informative.  
  # These are calcuated on R side and read in as data.

  # Informative priors on initial conditions based on first year's observation of family groups
  fg[1] ~ dpois(y[1])
  N[1] ~ dlnorm(log(y[1] / p),tau.p)

  # process model
  for (t in 2:(y.endyr + 1)) {  # the last year is a forecast with known harvest data
      mu[t] <- log(max(.0001, lambda * (N[t - 1] - y.H[t - 1])))
      N[t] ~ dlnorm(mu[t], tau.p)
      fg[t] <- N[t] * p
    }# end of process model
        
  # data model
  for (t in 2:y.endyr) {   
      y[t] ~ dpois(p * N[t])  
  }# end of data model

  # Model checking
  for (t in 1:y.endyr) {
    # simulate new data for posterior predicitve check
      y.rep[t] ~ dpois(p * N[t])
      # accumlate test statistics for posterior predictive check
      sq[t] <- (y[t] - p * N[t])^2
      sq.rep[t] <- (y.rep[t] - p * N[t])^2
    # compute residuals for autocorrelation check
    e[t] <- y[t] - fg[t]
  }

  # calculate Bayesian P value
  fit <- sum(sq[])
  fit.new <- sum(sq.rep[])
  pvalue <- step(fit.new - fit)

} #end of model

",fill=TRUE)
sink()
}



n.update = 10000
n.iter = 50000
n.adapt = 5000
n.thin = 1

jm = jags.model("LynxHarvest.R", data = data, inits = inits, n.adapt = n.adapt, n.chains = length(inits))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 18
##    Unobserved stochastic nodes: 43
##    Total graph size: 309
## 
## Initializing model
update(jm, n.iter = n.update)
z = coda.samples(jm, variable.names = c("lambda", "sigma.p", "p", "e", "N", "fg", "pvalue", "fit", "fit.new"), 
  n.iter = n.iter, thin = n.thin)


  1. Check convergence of MCMC chains for model parameters, process variance, and latent states. This will probably require using the excl option in MCMCsummary.


MCMCsummary(z, excl = c("fit.new", "fit", "e"), n.eff = TRUE, digits = 2)
##           mean     sd    2.5%    50%  97.5% Rhat  n.eff
## N[1]    500.00 77.000 370.000 490.00 670.00    1  13269
## N[2]    470.00 57.000 370.000 460.00 590.00    1   8992
## N[3]    390.00 45.000 310.000 390.00 490.00    1   9046
## N[4]    310.00 38.000 240.000 310.00 390.00    1   9196
## N[5]    270.00 34.000 210.000 270.00 350.00    1   9286
## N[6]    240.00 32.000 190.000 240.00 310.00    1   9163
## N[7]    230.00 31.000 170.000 230.00 290.00    1   9235
## N[8]    250.00 32.000 190.000 250.00 320.00    1   9620
## N[9]    270.00 33.000 210.000 260.00 340.00    1   9414
## N[10]   270.00 34.000 210.000 270.00 350.00    1   9693
## N[11]   300.00 37.000 240.000 300.00 380.00    1  10094
## N[12]   300.00 36.000 240.000 300.00 380.00    1   9601
## N[13]   240.00 30.000 190.000 240.00 300.00    1   9602
## N[14]   220.00 28.000 170.000 220.00 280.00    1   9297
## N[15]   160.00 24.000 120.000 160.00 220.00    1   9455
## N[16]   140.00 22.000  95.000 130.00 180.00    1   8943
## N[17]   150.00 22.000 100.000 140.00 190.00    1   9762
## N[18]   180.00 26.000 140.000 180.00 240.00    1  10955
## N[19]   210.00 34.000 150.000 200.00 280.00    1  12773
## N[20]   180.00 53.000 100.000 170.00 310.00    1  26622
## fg[1]    82.00  9.100  65.000  82.00 100.00    1 151076
## fg[2]    77.00  7.200  64.000  76.00  92.00    1  33337
## fg[3]    64.00  5.700  53.000  64.00  76.00    1  62144
## fg[4]    51.00  5.000  42.000  51.00  61.00    1  48742
## fg[5]    45.00  4.600  36.000  45.00  54.00    1  42819
## fg[6]    40.00  4.300  32.000  40.00  49.00    1  34891
## fg[7]    38.00  4.300  29.000  38.00  46.00    1  25997
## fg[8]    41.00  4.300  32.000  40.00  49.00    1  47993
## fg[9]    44.00  4.500  35.000  44.00  53.00    1  53007
## fg[10]   45.00  4.600  36.000  45.00  54.00    1  57045
## fg[11]   50.00  4.900  41.000  49.00  60.00    1  48086
## fg[12]   49.00  4.800  41.000  49.00  59.00    1  47418
## fg[13]   39.00  4.000  32.000  39.00  48.00    1  48971
## fg[14]   36.00  3.700  29.000  36.00  44.00    1  43200
## fg[15]   27.00  3.300  21.000  27.00  34.00    1  34470
## fg[16]   22.00  3.200  16.000  22.00  28.00    1  17709
## fg[17]   24.00  3.200  18.000  24.00  30.00    1  24900
## fg[18]   30.00  3.700  24.000  30.00  38.00    1  36599
## fg[19]   34.00  4.900  25.000  33.00  45.00    1  30391
## fg[20]   30.00  8.300  18.000  29.00  50.00    1  44007
## lambda    1.10  0.048   0.970   1.10   1.20    1  42437
## p         0.16  0.012   0.140   0.16   0.19    1   4862
## pvalue    0.58  0.490   0.000   1.00   1.00    1 100165
## sigma.p   0.17  0.061   0.077   0.16   0.31    1   9719


  1. Conduct posterior predictive checks by simulating a new dataset for family groups (\(f_t\)) at every MCMC iteration. Calculate a Bayesian p value using the sums of squared discrepancy between the observed and the predicted number of family groups based on observed and simulated data,

\[T^{observed} = \sum_{t=1}^{n}(f_{t}^{observed}-N_{t}\phi)^{2} \\T^{model} = \sum_{t=1}^{n}(f_{t}^{simulated}-N_{t}\phi)^{2}.\] The Bayesian p value is the proportion of MCMC iterations for which \(T^{model}>T^{obs}\).


par(mfrow = c(1, 1))
fit.new = MCMCchains(z, "fit.new")
fit = MCMCchains(z, "fit")
plot(fit.new,fit, xlab = "Discrepancy observed", ylab = "Discrepancy simulated", cex = .05, xlim = c(0,3000), 
  ylim = c(0,3000))
abline(0,1)
p = MCMCpstr(z, params = "pvalue", func = mean)
text(500, 2500, paste("Bayesian P value = ", as.character(signif(p$pvalue, 2))))


Assure yourself that the process model adequately accounts for temporal autocorrelation in the residuals— allowing the assumption that they are independent and identically distributed. To do this, include a derived quantity

\[e_t=y_t-N_t\phi,\]

in your JAGS code and JAGS object. Use the following code or something like it to examine how autocorrelation in the residuals changes with time lag.

  1. Write a paragraph describing how to interpret the plot produced by this function.
acf(unlist(MCMCpstr(z, param = "e", func = mean)), main = "", lwd = 3, ci = 0)


The autocorrelation of time series is the Pearson correlation between values of the process at two different times, as a function of the lag between the two times. It can take on values between -1 and 1. Values close to 1 or -1 indicate a high degree of correlation. The residuals are not auto correlated if their values drop close to 0 at relatively short lags and alternate between positive and negative values at subsequent lags. This plot reveals no autocorrelation in the residuals of our model.


  1. Plot the median of the marginal posterior distribution of the number of lynx family groups over time (1998-2016) including a highest posterior density interval. Include your forecast for 2017 (the predictive process distribution) in this plot.


# Plot true, unobserved state and forecast vs observations.
par(mfrow = c(1, 1))
# +1 gives the one year forecast
years = seq(1998, y[nrow(y), 1] + 1)
fg = MCMCpstr(z, params = "fg", func = median)
fgHPDI = MCMCpstr(z, params = "fg", func = function(x) hdi(x, .95))
y2 = c(y$census, NA)
plot(years,y2, ylim = c(0, 100), ylab = "Number of Lynx Family Groups", xlab = "Years")
lines(years, fg$fg)
lines(years, fgHPDI$fg[,2], lty = "dashed")
lines(years, fgHPDI$fg[,1], lty = "dashed")
abline(h=26, lty = "dotted", col = "red")
abline(h=32, lty = "dotted", col = "red")
Median population size of lynx (solid line) during 1997-2016 and forecasts for 2017 with 95% credible intervals (dashed lines). Red dotted lines give acceptable range of number of family groups determined in public input process.

Median population size of lynx (solid line) during 1997-2016 and forecasts for 2017 with 95% credible intervals (dashed lines). Red dotted lines give acceptable range of number of family groups determined in public input process.


  1. ADVANCED Lynx harvest decisions are made before the population is censused, even though harvest actually occurs shortly after the census. This mismatch occurs because of the time that it takes to properly issue hunting permits/licenses. Make a forecast of the number of family groups in 2018 assuming five alternative levels for 2017 harvest (0, 10, 25, 50, and 75 animals). Environmentalists and hunters have agreed on a acceptable range for lynx abundance in the unit, 26 - 32 family groups. Compute the probability that the post-harvest number of family groups will be below, within, and above this range during 2018. Tabulate these values. Hint: Set up a “model experiment” in your JAGS code where you forecast the number of lynx family groups during 2018 under the specified levels of harvest. Extract the MCMC chains for the forecasted family groups( e.g., fg.hat) using MCMCchains Use the ecdf function on the R side to compute the probabilities that the forecasted number groups will be below, within, or above the acceptable range.


# Levels of  Harvest to evaluate relative to goals for forecasting part.
h = c(0, 10, 25, 50, 75)

data = list(
    y.endyr = endyr,
    y.a = shapes[1], 
    y.b = shapes[2],
    y.H = y$harvest,
    y = y$census,
    h = h)

inits = list(
    list(lambda = 1.2, sigma.p = .01, N = n),
    list(lambda = 1.01,sigma.p = .2, N = n * 1.2),
    list(lambda = .95, sigma.p = .5, N = n * .5))
{
sink("LynxHarvest.R")
cat("

model{

  # priors
  sigma.p ~ dunif(0, 5)
  tau.p <- 1 / sigma.p^2
  lambda ~ dunif(0, 10)
  p ~ dbeta(y.a, y.b)  

  # Get parameters y.a and y.b from mean and sd using moment matching to make this prior informative.  
  # These are calcuated on R side and read in as data.

  # Informative priors on initial conditions based on first year's observation of family groups
  fg[1] ~ dpois(y[1])
  N[1] ~ dlnorm(log(y[1] / p),tau.p)

  # process model
  for (t in 2:(y.endyr + 1)) {  # the last year is a forecast with known harvest data
      mu[t] <- log(max(.0001, lambda * (N[t - 1] - y.H[t - 1])))
      N[t] ~ dlnorm(mu[t], tau.p)
      fg[t] <- N[t] * p
    }# end of process model
        
  # data model
  for (t in 2:y.endyr) {   
      y[t] ~ dpois(p * N[t])  
  }# end of data model

  # Model checking
  for (t in 1:y.endyr) {
    # simulate new data for posterior predicitve check
      y.rep[t] ~ dpois(p * N[t])
      # accumlate test statistics for posterior predictive check
      sq[t] <- (y[t] - p * N[t])^2
      sq.rep[t] <- (y.rep[t] - p * N[t])^2
    # compute residuals for autocorrelation check
    e[t] <- y[t] - fg[t]
  }

  # calculate Bayesian P value
  fit <- sum(sq[])
  fit.new <- sum(sq.rep[])
  pvalue <- step(fit.new - fit)

  # Forecast effects of different harvest regeimes on number of family grops during 2018
    for(i in 1:length(h)) {
        mu.hat[i] <- log(max(.001, lambda * (N[y.endyr + 1] - h[i]))) 
        N.hat[i] ~ dlnorm(mu.hat[i], tau.p) # Nhat forecasts 2 years out
        fg.hat[i] <- N.hat[i] * p
  }
    
} #end of model

",fill=TRUE)
sink()
}
n.update = 10000
n.iter = 50000
n.adapt = 5000
n.thin = 1

jm = jags.model("LynxHarvest.R", data = data, inits = inits, n.adapt = n.adapt, n.chains = length(inits))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 18
##    Unobserved stochastic nodes: 48
##    Total graph size: 345
## 
## Initializing model
update(jm, n.iter = n.update)
z2 = coda.samples(jm, variable.names = "fg.hat", n.iter = n.iter, thin = n.thin)
# Acceptable limits on poplation size, determined by public input process.
lower = 26
upper = 32
p.in = numeric(length(h))
p.over = numeric(length(h))
p.under = numeric(length(h))

# get chains for harvest experiment
fg.hat = MCMCchains(z2, "fg.hat")

for (j in 1:length(h)) {
    p1 = ecdf(fg.hat[, j])(upper)
    p.under[j] = ecdf(fg.hat[, j])(lower)
    p.in[j] = p1 - p.under[j]
    p.over[j] = 1 - p1
}

# trim to reasonable signficiant digits
p.under = signif(p.under, 2)
p.in = signif(p.in, 2)
p.over = signif(p.over, 2)

alt.table = as.data.frame(cbind(h, p.under, p.in, p.over))
names(alt.table) = c("Harvest", "P(under)", "P(in)", "P(over)")
alt.table
##   Harvest P(under) P(in) P(over)
## 1       0     0.31 0.270   0.430
## 2      10     0.39 0.250   0.360
## 3      25     0.51 0.220   0.270
## 4      50     0.69 0.150   0.160
## 5      75     0.81 0.094   0.095