# Point-wise log-likelihood calculation for CJS with logistic predictor

Hiya folks! First post, so apologies if this isn’t articulated well. I’m working with a Cormack-Jolly-Seber model using the marginalization method outlined in Yackulic et al. (2020). I’ve adapted it to include a predictor of the survival rates (e.g. some climatic variable). Additionally, I’ve added some code to set the detection probablity of unsampled years to 0. However, I’m uncertain of how to extract the point-wise log-likelihood values for use in loo. I suppose I need some vector in a generated quantity block that is the likelihood of each capture history given my estimates for the survival and detection probabilities. Then, I think I would append this vector with the likelihood of these survival estimates given the values of my predictor variable, but I am less sure about this. Overall, I’m a bit stumped and am wondering if others have done anything similar or have ideas. My code without this generated quantity is below.

``````// CJS Model with time varying survival (s) and capture probability (p)

data{
int<lower = 1> NsumCH; // number of capture histories
int<lower = 1> n_occasions; // number of capture occasions
int<lower = 1, upper = 2> sumCH[NsumCH, n_occasions]; // capture histories
int<lower = 1> sumf[NsumCH]; // first capture occasion
int<lower = 1> sumFR[NsumCH]; // frequency of each capture history
int<lower = 1> n_missing; // number of unsampled years
int<lower = 1> missing[n_missing]; // unsampled years
int<lower = 1> n_observed; // number of sampled years
int<lower = 1> observed[n_observed]; // sampled years
vector[n_occasions-1] pred; // predictor variable each year
}

parameters{
real<lower = 0, upper = 1> free_p[n_occasions - 1]; // detection probability for all years
real<lower = 0, upper = 1> logit_phi[n_occasions - 1]; // logit survival probability
real phi_0; // intercept for survival model
real phi_pred; // effect of  pred on survival
}

transformed parameters{
simplex tr[2,n_occasions - 1]; // survival likelihoods for marginalization
simplex pmat[2,n_occasions - 1]; // detection likelihoods for marginalization
// real mu[n_occasions-1];
real<lower = 0, upper = 1> p[n_occasions - 1]; // detection probability with unsampled years = 0
real<lower=0,upper=1> phi[n_occasions-1]; // survival probability

for(i in missing){
p[i] = 0; // set unsampled p to 0
}

for(i in observed){
p[i] = free_p[i]; // fill in sampled p with a real estimate
}

for(i in 1:n_occasions - 1){
phi[i] = inv_logit(phi_0+pred[i]*phi_pred); // survival is a function of predictor vars
}

for(k in 1:n_occasions - 1){
tr[1,k,1] = phi[k];
tr[1,k,2] = (1 - phi[k]);
tr[2,k,1] = 0;
tr[2,k,2] = 1;

pmat[1,k,1] = p[k];
pmat[1,k,2] = (1 - p[k]);
pmat[2,k,1] = 0;
pmat[2,k,2] = 1;
}
}

model{
vector pz[n_occasions]; // marginalized likelihoods

free_p ~ uniform(0,1);

for(i in 1:NsumCH){
pz[sumf[i],1] = 1;
pz[sumf[i],2] = 0;

for(k in (sumf[i] + 1):n_occasions){
pz[k,1] = pz[k-1,1] * tr[1,k-1,1] * pmat[1,k-1,sumCH[i,(k)]];
pz[k,2] = (pz[k-1, 1] * tr[1,k-1,2] + pz[k-1, 2]) * pmat[2,k-1,sumCH[i,(k)]];
}

target += sumFR[i] * log(sum(pz[n_occasions]));

}

// for(i in 1:n_occasions-1){
//   target += normal_lpdf(logit_phi[i] | mu[i], 10);
// }
}
``````

Any help would be great! Let me know if I’m missing important information in this post or if I’m demonstrating a misconception or gap in my knowledge. Thanks in advance!

Hi and welcome! Unfortunately, I don’t have time this week to help with this question in detail, but you might find this writing helpful: Occupancy model likelihoods: advanced topics

It’s written about occupancy models, not CJS models, but it elaborates on some concepts about the pointwise likelihood that transfer over and might help out.

That link was very helpful! If I understand correctly, and am applying your occupancy model example correctly to my CJS model, I need to calculate the likelihood of each observation (I assume each observation after the initial observation) per animal, given the the survival and detection probabilities. So with my code, I think it would look like:

``````
generated quantities {
vector[n_lik] log_lik;
int counter;
counter = 1;
for(i in 1:NsumCH){ // doesn't include critters observed first on last occasion
if (sumf[i] >= n_occasions){
}
else{
for(j in sumf[i]+1:n_occasions){
for(k in 1:sumFR[i]){
log_lik[counter] = bernoulli_lpmf(CH[i,j-1]|phi[j-1]*p[j-1]);
counter += 1;
}
}
}
}
}
``````

with the n_lik being the expected number of pointwise likelihoods I STAN should be returning.