# Generated quantities and divergent chains

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
}

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 = beta1;
beta = beta2;
beta = 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) {
p[i, t] = inv_logit(beta[group[i]] + gamma[t]*p_adjust[t]);

}
}

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 + gamma .* p_adjust); // Back-transformed detection of tag delay groups
p_g3 = inv_logit(beta + 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
}

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 = beta1;
beta = beta2;
beta = 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) {
p[i, t] = inv_logit(beta[group[i]] + gamma[t]);

}
}

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 + gamma); // Back-transformed detection of tag delay groups
p_g3 = inv_logit(beta + 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]]);
}
}
}

``````

I don’t think the two models are the same.

The first is:

``````p[i, t] = inv_logit(beta[group[i]] + gamma[t]*p_adjust[t]);
``````

The second is:

``````p[i, t] = inv_logit(beta[group[i]] + gamma[t]);
...
y[i, t] ~ bernoulli(p[i, t - 1] * p_adjust[t - 1]);
``````

For them to be the same the first would need to be:

``````p[i, t] = inv_logit(beta[group[i]] + gamma[t]) * p_adjust[t];
``````

Regarding divergences, you should use `bernoulli_logit` (that includes the `inv_logit` for you) instead of `bernoulli` and `inv_logit` separately. That can help with numerics.

If you need `p[i, t]` on the non-logit scale, you can compute that separately in generated quantities.

1 Like

Thanks very much, things are looking a lot more sensible now!

1 Like