Hello again,
I recently had a minor syntax error that was kindly corrected for me see here. The code will now run, but I get a lot of divergent chains and weird estimates.
Basically I have a simple mark-recapture model and want to back-transform some estimates of time-varying detection per group, but because I have to adjust for uneven detection effort I have an additional term p_adjust that I need to use to correct estimates during periods I know detection was less than 1 or in some cases 0. The thing that is confusing me is that if I take the p_adjust out of the model constraints and out of the generated quantities block and just include it directly in the likelihood (the second model pasted below) the model appears to run fine. Anyone able to enlighten me as to what is going here? I thought, obviously wrongly given my lack of understanding, that they should be vaguely comparable.
Any insight would be much appreciated.
Many thanks.
First model that results in divergent chains…
functions {
int first_capture(int[] y_i) {
for (k in 1:size(y_i))
if (y_i[k])
return k;
return 0;
}
int last_capture(int[] y_i) {
for (k_rev in 0:(size(y_i) - 1)) {
int k = size(y_i) - k_rev;
if (y_i[k])
return k;
}
return 0;
}
matrix prob_uncaptured(int nind, int n_occasions,
matrix p, matrix phi) {
matrix[nind, n_occasions] chi;
for (i in 1:nind) {
chi[i, n_occasions] = 1.0;
for (t in 1:(n_occasions - 1)) {
// Compoud declaration was enabled in Stan 2.13
int t_curr = n_occasions - t;
int t_next = t_curr + 1;
chi[i, t_curr] = (1 - phi[i, t_curr]) + phi[i, t_curr] * (1 - p[i, t_next - 1]) * chi[i, t_next];
}
}
return chi;
}
}
data {
int<lower=0> nind; // Number of individuals
int<lower=2> n_occasions; // Number of capture occasions
int<lower=0,upper=1> y[nind, n_occasions]; // Capture-history
int<lower=1> g; // Number of groups
int<lower=1,upper=g> group[nind]; // Groups
real<lower=0,upper=1>p_adjust[n_occasions-1];
int<lower=2, upper=4> phi_adjust[n_occasions-1];
}
transformed data {
int n_occ_minus_1 = n_occasions - 1;
int<lower=0,upper=n_occasions> first[nind];
int<lower=0,upper=n_occasions> last[nind];
real beta1 = 0; // Corner constraint
for (i in 1:nind)
first[i] = first_capture(y[i]);
for (i in 1:nind)
last[i] = last_capture(y[i]);
}
parameters {
real<lower=0,upper=1> mean_phi; // Mean survival
real<lower=0,upper=1> mean_p; // Mean recapture
vector[n_occ_minus_1] gamma; // Time effects
//vector<lower=0,upper=1>[g] p_g; // Group-spec. recapture
real beta2; // Prior for difference in tag delay
real beta3; // Prior for difference in tag delay
}
transformed parameters {
matrix<lower=0,upper=1>[nind, n_occ_minus_1] phi;
matrix<lower=0,upper=1>[nind, n_occ_minus_1] p;
matrix<lower=0,upper=1>[nind, n_occasions] chi;
vector[g] beta;
beta[1] = beta1;
beta[2] = beta2;
beta[3] = beta3;
// Constraints
for (i in 1:nind) {
for (t in 1:(first[i] - 1)) {
phi[i, t] = 0;
p[i, t] = 0;
}
for (t in first[i]:n_occ_minus_1) {
phi[i, t] = exp(phi_adjust[t]*(log(mean_phi)/3));
// phi_adjust adjusts for uneven sampling periods in April 2012
p[i, t] = inv_logit(beta[group[i]] + gamma[t]*p_adjust[t]);
//p_adjust adjusts for uneven detection (recapture) effort
}
}
chi = prob_uncaptured(nind, n_occasions, p, phi);
}
model {
// Priors
// Uniform priors are implicitly defined.
// mean_phi ~ uniform(0, 1);
// mean_p ~ uniform(0, 1);
// p_g ~ uniform(0, 1);
beta2 ~ normal(0, 10)T[-10,10];
beta3 ~ normal(0, 10)T[-10,10];
gamma ~ normal(0, 10);
// Likelihood
for (i in 1:nind) {
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
1 ~ bernoulli(phi[i, t - 1]);
y[i, t] ~ bernoulli(p[i, t - 1]);
}
1 ~ bernoulli(chi[i, last[i]]);
}
}
}
generated quantities {
vector<lower=0,upper=1>[n_occ_minus_1] p_g1;
vector<lower=0,upper=1>[n_occ_minus_1] p_g2;
vector<lower=0,upper=1>[n_occ_minus_1] p_g3;
vector[nind] log_lik;
// This is where I'm getting the error vector*real array
p_g1 = inv_logit(gamma .* p_adjust); // Back-transformed detection of tag delay groups
p_g2 = inv_logit(beta[2] + gamma .* p_adjust); // Back-transformed detection of tag delay groups
p_g3 = inv_logit(beta[3] + gamma .* p_adjust); // Back-transformed detection of tag delay groups
for (i in 1:nind) {
log_lik[i] = 0;
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
log_lik[i] = log_lik[i] + bernoulli_lpmf(1|phi[i, t - 1]);
log_lik[i] = log_lik[i] + bernoulli_lpmf(y[i, t]|p[i, t - 1]);
}
log_lik[i] = log_lik[i] + bernoulli_lpmf(1|chi[i, last[i]]);
}
}
}
second model where p_adjust is included directly in the likelihood
functions {
int first_capture(int[] y_i) {
for (k in 1:size(y_i))
if (y_i[k])
return k;
return 0;
}
int last_capture(int[] y_i) {
for (k_rev in 0:(size(y_i) - 1)) {
int k = size(y_i) - k_rev;
if (y_i[k])
return k;
}
return 0;
}
matrix prob_uncaptured(int nind, int n_occasions,
matrix p, matrix phi) {
matrix[nind, n_occasions] chi;
for (i in 1:nind) {
chi[i, n_occasions] = 1.0;
for (t in 1:(n_occasions - 1)) {
// Compoud declaration was enabled in Stan 2.13
int t_curr = n_occasions - t;
int t_next = t_curr + 1;
chi[i, t_curr] = (1 - phi[i, t_curr]) + phi[i, t_curr] * (1 - p[i, t_next - 1]) * chi[i, t_next];
}
}
return chi;
}
}
data {
int<lower=0> nind; // Number of individuals
int<lower=2> n_occasions; // Number of capture occasions
int<lower=0,upper=1> y[nind, n_occasions]; // Capture-history
int<lower=1> g; // Number of groups
int<lower=1,upper=g> group[nind]; // Groups
real<lower=0,upper=1>p_adjust[n_occasions-1];
int<lower=2, upper=4> phi_adjust[n_occasions-1];
}
transformed data {
int n_occ_minus_1 = n_occasions - 1;
int<lower=0,upper=n_occasions> first[nind];
int<lower=0,upper=n_occasions> last[nind];
real beta1 = 0; // Corner constraint
for (i in 1:nind)
first[i] = first_capture(y[i]);
for (i in 1:nind)
last[i] = last_capture(y[i]);
}
parameters {
real<lower=0,upper=1> mean_phi; // Mean survival
real<lower=0,upper=1> mean_p; // Mean recapture
vector[n_occ_minus_1] gamma; // Time effects
//vector<lower=0,upper=1>[g] p_g; // Group-spec. recapture
real beta2; // Prior for difference in tag delay
real beta3; // Prior for difference in tag delay
}
transformed parameters {
matrix<lower=0,upper=1>[nind, n_occ_minus_1] phi;
matrix<lower=0,upper=1>[nind, n_occ_minus_1] p;
matrix<lower=0,upper=1>[nind, n_occasions] chi;
vector[g] beta;
beta[1] = beta1;
beta[2] = beta2;
beta[3] = beta3;
// Constraints
for (i in 1:nind) {
for (t in 1:(first[i] - 1)) {
phi[i, t] = 0;
p[i, t] = 0;
}
for (t in first[i]:n_occ_minus_1) {
phi[i, t] = exp(phi_adjust[t]*(log(mean_phi)/3));
// phi_adjust adjusts for uneven sampling periods in April 2012
p[i, t] = inv_logit(beta[group[i]] + gamma[t]);
//p_adjust adjusts for uneven detection (recapture) effort
}
}
chi = prob_uncaptured(nind, n_occasions, p, phi);
}
model {
// Priors
// Uniform priors are implicitly defined.
// mean_phi ~ uniform(0, 1);
// mean_p ~ uniform(0, 1);
// p_g ~ uniform(0, 1);
beta2 ~ normal(0, 10)T[-10,10];
beta3 ~ normal(0, 10)T[-10,10];
gamma ~ normal(0, 10);
// Likelihood
for (i in 1:nind) {
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
1 ~ bernoulli(phi[i, t - 1]);
y[i, t] ~ bernoulli(p[i, t - 1] * p_adjust[t -1]);
}
1 ~ bernoulli(chi[i, last[i]]);
}
}
}
generated quantities {
vector<lower=0,upper=1>[n_occ_minus_1] p_g1;
vector<lower=0,upper=1>[n_occ_minus_1] p_g2;
vector<lower=0,upper=1>[n_occ_minus_1] p_g3;
vector[nind] log_lik;
// This is where I'm getting the error vector*real array
p_g1 = inv_logit(gamma); // Back-transformed detection of tag delay groups
p_g2 = inv_logit(beta[2] + gamma); // Back-transformed detection of tag delay groups
p_g3 = inv_logit(beta[3] + gamma); // Back-transformed detection of tag delay groups
for (i in 1:nind) {
log_lik[i] = 0;
if (first[i] > 0) {
for (t in (first[i] + 1):last[i]) {
log_lik[i] = log_lik[i] + bernoulli_lpmf(1|phi[i, t - 1]);
log_lik[i] = log_lik[i] + bernoulli_lpmf(y[i, t]|p[i, t - 1] * p_adjust[t - 1]);
}
log_lik[i] = log_lik[i] + bernoulli_lpmf(1|chi[i, last[i]]);
}
}
}