Firstly apologies for writing a very similar post on something that has already come up a number of times.
I’m using a mark-recapture model, where I have to adjust for periods of uneven recapture effort, this information is included in the data block as p_adjust, and subsequently included in the model constraints. I then need to back-transform the estimates in the generated quantities block, and this is where I get the error.
As the model stands I’m trying to multiply a vector by a real array which is not allowed, but neither is multiplying a vector by a vector, so I know what my problem is but not the best way to correct it. I’m sure this is a relatively straight forward issue, and I’d be very grateful if someone more knowledge could suggest the easiest way to define p_adjust that will avoid this error.
Many thanks!
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(gamma + beta[2]*p_adjust); // Back-transformed detection of tag delay groups
p_g3 = inv_logit(gamma + beta[3]*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]]);
}
}
}