In cttedwards/bdm: Bayesian biomass dynamics model

knitr::opts_chunk$set(fig.path = 'fig/bdm-newmodel-', fig.width=  6, tidy = FALSE, message = FALSE, warning = FALSE, collapse = TRUE, comment = "#>")

Example applications of non-default (i.e. user-specified) models in the bdm R-package are given here, based on data from fisheries in New Zealand.

library(bdm)

Chatham Rise Hake

The model is fitted to data from the chatham rise hake fishery in New Zealand, which consists of catches, a commerical abundance index and a survey index. The data are used to initialise an empirical data object (bdmData) that, by default, renormalises the indices to a geometric mean of one. See ?bdmData for details.

data(haknz)
dat <- bdmData(harvest = haknz$landings,index = cbind(haknz$survey, haknz$cpue), time = haknz$year, renormalise = TRUE)
plot(dat)

The revised model could be coded with the following changes relative to the default model:

  new_model <- "
    ...
    parameters {
      real<lower=3,upper=30> logK;
      real<lower=0,upper=2> r;
      real<lower=0> x[T];
      real<lower=0> sigmap2; 
    }
    ...
    model {

      // prior densities for
      // estimated parameters
      // ********************
      logK ~ uniform(3.0,30.0);
      r ~ lognormal(-1.0,0.20);
      sigmap2 ~ inv_gamma(0.001,0.001);

      ...
    }
    ...
  "

  new_model <- "
    data {
      int T;
      int I;
      real index[T,I];
      real harvest[T];
      real n;
      real sigmao[T,I];
      real sigmap;
    }
    parameters {
      real<lower=3,upper=30> logK;
      real<lower=0,upper=2> r;
      real<lower=0> x[T];
      real<lower=0> sigmap2; 
    }
    transformed parameters {

      real q[I];

      // variance terms
      real sigmao2[T,I];

      // fletcher-schaefer
      // parameters
      real dmsy;
      real h;
      real m;
      real g;

      dmsy <- pow((1/n),(1/(n-1)));
      h <- 2*dmsy;
      m <- r*h/4;
      g <- pow(n,(n/(n-1)))/(n-1);

      // variance terms
      for(t in 1:T)
        for(i in 1:I)
            sigmao2[t,i] <- square(sigmao[t,i]);

      // compute mpd catchability assuming 
      // variable sigmao over time assuming
      // uniform prior on q
      {
        real sum1;
        real sum2;
        real p;
        for(i in 1:I){
          sum1 <- 0.0;
          sum2 <- 0.0;
          p <- 0.0;
          for(t in 1:T){
            if(index[t,i]>0.0 && x[t]>0.0) {
              sum1 <- sum1 + log(index[t,i]/x[t])/sigmao2[t,i];
              sum2 <- sum2 + 1/sigmao2[t,i];
              p <- p + 1.0;
            }
          }
          if(p>2.0) { q[i] <- exp((0.5 * p + sum1) / sum2);
            } else q[i] <- 0.0;
        }
      }
    } 
    model {

      // prior densities for
      // estimated parameters
      // ********************
      logK ~ uniform(3.0,30.0);
      r ~ lognormal(-1.0,0.20);
      sigmap2 ~ inv_gamma(0.001,0.001);

      // state equation
      // **************
      {
        real H;
        real mu;
        x[1] ~ lognormal(log(1.0)-sigmap2/2, sqrt(sigmap2));
        for(t in 2:T) {
          H <- fmin(exp(log(harvest[t-1]) - logK),x[t-1]);
          if(x[t-1]<=dmsy) mu <- x[t-1] + r * x[t-1] * (1 - x[t-1]/h) - H;
          if(x[t-1]> dmsy) mu <- x[t-1] + g * m * x[t-1] * (1 - pow(x[t-1],(n-1))) - H;
          if(mu > 0.0) mu <- log(mu) - sigmap2/2;
          else mu <- log(0.01) - sigmap2/2;
          x[t] ~ lognormal(mu, sqrt(sigmap2));
        }
      }

      // observation equation
      // ********************
      {
        real mu;
        for(i in 1:I){
          for(t in 1:T){
            if(index[t,i]>0.0 && x[t]>0.0 && q[i]>0.0) {
              mu <- log(q[i]*x[t]) - sigmao2[t,i]/2;
              index[t,i] ~ lognormal(mu,sigmao[t,i]);
            }
          }
        }
      }

      // apply penalty for H>0.95
      // ************************
      for(t in 1:T){
      real H_; H_ <- harvest[t]/exp(log(x[t]) + logK);
        if(H_>0.95) {
          increment_log_prob(-log(H_/0.95) * (1/sigmap2));
        }
      }
    }
    generated quantities {

      real biomass[T];
      real depletion[T];
      real harvest_rate[T];
      real surplus_production[T];

      real epsilon_o[T,I];
      real epsilon_p[T];

      real current_biomass;
      real current_depletion;
      real current_harvest_rate;

      real msy;
      real biomass_at_msy;
      real harvest_rate_at_msy;

      real observed_index[T,I];
      real predicted_index[T,I];

      {
        real H;
        for(t in 2:T) {
          H <- fmin(exp(log(harvest[t-1]) - logK),x[t-1]);
          if(x[t-1]<=dmsy) epsilon_p[t-1] <- x[t]/(x[t-1] + r * x[t-1] * (1 - x[t-1]/h) - H);
          if(x[t-1]> dmsy) epsilon_p[t-1] <- x[t]/(x[t-1] + g * m * x[t-1] * (1 - pow(x[t-1],(n-1))) - H);
        }
        epsilon_p[T] <- lognormal_rng(log(1.0)-sigmap2/2, sqrt(sigmap2));
      }

      for(t in 1:T) {
        biomass[t] <- x[t] * exp(logK);
        depletion[t] <- x[t];
        harvest_rate[t] <- harvest[t]/exp(log(x[t]) + logK);
        if(x[t]<=dmsy) surplus_production[t] <- r * x[t] * (1 - x[t]/h) * epsilon_p[t];
        if(x[t]> dmsy) surplus_production[t] <- g * m * x[t] * (1 - pow(x[t],(n-1))) * epsilon_p[t];
      }

      current_biomass <- biomass[T];
      current_depletion <- x[T];
      current_harvest_rate <- harvest_rate[T];

      msy <- m * exp(logK);
      biomass_at_msy <- dmsy * exp(logK);
      harvest_rate_at_msy <- m / dmsy;

      for(i in 1:I){
        for(t in 1:T){
          observed_index[t,i] <- index[t,i];
          predicted_index[t,i] <- q[i]*x[t];
        }
      }

      for(t in 1:T){
        for(i in 1:I){
          epsilon_o[t,i] <- observed_index[t,i]/predicted_index[t,i];
        }
      }
    }
  "

This new formulation includes a revision of our assumptions regarding the process error, allowing the variance to be estimated. We can update the model code using a list argument supplied to the updatePrior function. It is then necessary to compile the model:

mdl <- bdm(model.code = new_model)
mdl <- updatePrior(mdl, prior = list(par = 'r', meanlog = -0.91, sdlog = 0.40))
mdl <- compiler(mdl)

If a new model is loaded into the bdm object, and if the estimated parameters are not the same as the default, then the initialisation function within sampler() will not work. This function provides sensible values for the MCMC algorithm. If the new model were to be run without an initialisation function then initial parameter values will be generated randomly by rstan.

To construct a new inititialisation function we can use some of the auxilliary functions that are implemented for the default model:

init.r    <- getr(mdl)[['E[r]']]
init.logK <- getlogK(dat, r = init.r)
init.x    <- getx(dat, r = init.r, logK = init.logK)

These values can be used to create the initialisation function that must return a list of initial values:

init.func <- function() list(r = init.r, logK = init.logK, x = init.x, sigmap2 = 0.0025)
mdl <- sampler(mdl, dat, init = init.func)

histplot(mdl,par = c('r','logK','sigmap2'))

We note that the model is unable to resolve a value for logK, and estiamates an excessively high value for sigmap2. It would therefore be considered to have failed as an assessment. This can also be seen from the catchability traces, which behave poorly: