Pointwise log-likelihood for Mark-Recapture Model

Modeling
, ,
1

Hi,
I am trying to implement the individual Cormack-Jolly-Seber model provided in the Stan manual section 7.3 with different covariates and random effects. I would like to calculate LOO for each model with different plausible covariates for survival for model selection. I am having trouble trying to calculate the pointwise log-likelihood because of the indexing is based on the first and last capture in the matrix. See model block:

model {
  for (i in 1:I) {
    if (first[i] > 0) {
      for (t in (first[i]+1):last[i]) {
        1 ~ bernoulli(phi[t-1]);
        y[i, t] ~ bernoulli(p[t]);
      }
      1 ~ bernoulli(chi[last[i]]);
    }
  }
}

There is probably an easy trick to this but it is totally escaping me. I have tried setting individual log_lik components inside the loop and bring them together using append_row but I can’t seem to figure out how to increment the correct counter. For example, this is not working (note I’m omitting the variable declaration in generated quantities block):

generated quantities{
  for (i in 1:nind) {
    if (first[i] > 0) {
      for (t in (first[i] + 1):last[i]) {
        log_lik_1[(last[i] - (first[i] + 1) + 1) * (i-1) + (t - first[i])] = bernoulli_lpmf(1|phi[i, t - 1]);
        log_lik_2[(last[i] - (first[i] + 1) + 1) * (i-1) + (t - first[i])] = bernoulli_lpmf(y[i, t]|p[i, t - 1]);
      }
      log_lik_3[i] = bernoulli_lpmf(1|chi[i, last[i]]);
    }
  }
   log_lik_4 = append_row(log_lik_1, log_lik_2);
   log_lik = append_row(log_lik_3, log_lik_4);
}
}

The problem with the code above is that the for loop index with t is not triggered if the individual fish was only sampled once throughout the time series. Thus incrementing the log_lik_1 and log_lik_2 using the outer look for i is an incorrect index.

Any help or suggestions will be much appreciated.
Thanks,
Jason

2

I think I found a solution after some digging on this page. Apologies for not finding it before posting a duplicate question. The answer was supplied by user mbjoseph, thanks!

generated quantities{  
vector[nind] log_lik;
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]]);
    }
  }
}

Thanks,
Jason

2 Likes