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
  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]]);
      } 
    } 
}

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