Hello,
Tried to use the ecology tag (suggested here…https://stanecology.github.io/), but can’t find it. Anywhos:
I’m attempting to fit something like a multi-scale occupancy model (sensu Nichols et al. 2008). As the topic states, the problem I’ve faced is that that every iteration ends up with a divergent transition, and I’m trying to figure out whether this is fixable or not. I suspect the problem relates to a coding error/mis-specification of the likelihood, although it could be that the data in hand don’t support the model. At least, I have not found that increasing adapt.delta to 0.95 and upping the length of the warm-up to 3500 makes any difference.
Briefly, the complete data likelihood is like:
z_i~Bernoulli(Psi), where z_i denotes whether a site is occupied by some species
a_i_t~Bernoulli(theta), where a_i_t denotes whether the species is ‘available’ to be detected (e.g., it is in site i during time period t given that the site is ever occupied)
y_i_t_k~Bernoulli(z_ia_i_tp), where y_i_t_k is observed presence (with probability p conditional on occupancy and availability) or absence at one of multiple points k contained within site i at time t. The replicated observation points k are required for parameters p and theta to be identifiable.
The likelihood I’m employing in the stan program looks like so:
for (o in 1:nobs) {
//data is in long format, with a corresponding vector denoting site, time, and k indexing
if (obstate[o]==1) {
//known occupancy, known availability--i.e., species observed at site i at time t at at least 1 k
1~bernoulli(psi[site_idx[o]]);
1~bernoulli(theta[site_idx[o], Day[o]]);
y[o]~bernoulli(p[k_idx[o]]);
}
else if (obstate[o]==2){
//known occupancy, because observed at site i at least one time.
//unknown availability. Could have been unavailable, or available and not detected.
1~bernoulli(psi[site_idx[o]]);
target +=log_sum_exp(bernoulli_lpmf(1|theta[site_idx[o], Day[o]])+
bernoulli_lpmf(0 |p[k_idx[o]]),
bernoulli_lpmf(0|theta[site_idx[o], Day[o]]));
}
else {
//unknown occupancy and availability--never seen at site i.
//prob not occupied, prob occupied but not available, prob occupied & available but not detected
vector[3] lp;
lp[1]=bernoulli_lpmf(0|psi[site_idx[o]]);
lp[2]=bernoulli_lpmf(1|psi[site_idx[o]])+bernoulli_lpmf(0|theta[site_idx[o], Day[o]]);
lp[3]=bernoulli_lpmf(1|psi[site_idx[o]])+bernoulli_lpmf(1|theta[site_idx[o], Day[o]])+
bernoulli_lpmf(0|p[k_idx[o]]);
target+=log_sum_exp(lp);
}
} //end o loop
I’ve looked at this a lot, and have not seen anything obviously incorrect (could be more efficient), but maybe there’s an error in there.
The other thing that I thought might be an issue relates to some of the random effects structures employed. I think I’ve coded an uncentered model correctly, but perhaps I’m off? e.g.:
transformed parameters {
//other junk....
vector[365] a0_theta;
a0_theta[1]=a0_theta_raw[1];
for (t in 2:365){
a0_theta[t]=a0_theta[t-1] + a0_theta_raw[t]*var_a0_theta;
}
}//end transformed
model{
var_a0_theta ~ normal (0, 2.5); //<lower=0>
a0_theta_raw[1] ~ normal (0, 2.5);
for (t in 2:365){
a0_theta_raw[t] ~ std_normal();
}
}
****
When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use
```stan
model {
vector[N] mu = alpha + beta * x;
y ~ normal(mu, sigma);
}
instead of
model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}
To include mathematical notation in your post put LaTeX syntax between two $
symbols, e.g.,
p(\theta | y) \propto p(\theta) p(y | \theta).