I’m fitting about as simple of a model as there comes in population dynamics. It’s a logistic population model with observation and process error, where we have observed catches and an index of abundance from a fishery, and want to estimate population parameters such that given the catch history we can fit the abundance index. But I’m having some trouble making it converge properly across a broad range of fisheries. Specifically, while I can get a good fit (and convergence) for one fishery with some messing with the hyperparameters, when I try and use said model to fit data from many different fisheries the thing starts to run into all kinds of divergence and max_treedepth problems. The answer may simply be that I need to develop some routines for adjusting the hyperparameters for specific fisheries, or just be more patient.
I’m curious though if the community has any recommendations on ways to improve the overall structure of the model. In particular, I’ve never actually come across a dynamic population model on any of the Stan forums, and so am wondering if there are some tricks/nuances to using Stan in these circumstances (i.e. fitting a population model that evolves over time) that I should be using. I’ve included a link to the data and the stanfit object. The observation error (sigma_obs
) appears to be a particular problem and I’m wondering if there’s a better transformation here.
https://drive.google.com/open?id=1A4a185kDUMHnu3tYeqkAF_4eCIyA2-l4
data {
int<lower = 1> nt; // number of time steps to run
vector[nt] catch_t; // observed catch
vector[nt] log_index; // observed log population index
}
parameters {
real log_r; // log intrinsic growth rate
real<lower = 1> k_mult; // multiplier on mean catch to define carrying capacity
real<upper = 0> log_q; // catchability coefficient
real<lower = 0> sigma_obs; // observation error
real<lower = 0> sigma_r; // process error
vector[nt] log_rec_dev; // recruitment deviates
}
transformed parameters{
real k;
real r;
real q;
real temp;
real penalty;
vector[nt] b_t;
vector[nt] index_t;
vector[nt] rec_dev;
rec_dev = exp(sigma_r * log_rec_dev - sigma_r^2/2);
q = exp(log_q);
r = exp(log_r);
k = k_mult * (mean(catch_t) * 4 / r);
b_t = rep_vector(1e-3, nt);
index_t = rep_vector(1e-3, nt);
b_t[1] = k * rec_dev[1];
penalty = 1e-3;
for (t in 2:nt){
temp = ((b_t[t - 1] + b_t[t - 1]*r*(1 - b_t[t - 1] / k)) - catch_t[t - 1]) * rec_dev[t];
if (temp <= 0){ // if catches produce negative population
penalty += temp^2;
b_t[t] = 1e-3;
} else {
b_t[t] = temp;
}
} // close time loop
index_t = b_t * q; // calculate estimated population index
}
model {
log_index ~ normal(log(index_t), sigma_obs);
log_r ~ normal(log(0.2),2);
k_mult ~ normal(2,5);
log_q ~ normal(-6,0.5);
sigma_obs ~ cauchy(0,2.5);
sigma_r ~ normal(0.2,0.001); // constraining tightly for now
log_rec_dev ~ normal(0,1);
target += -penalty; // penalty for crashing the population
}