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

Modeling
1

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[2] tr[2,n_occasions - 1]; // survival likelihoods for marginalization
  simplex[2] 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[2] 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!

2

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.

3

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.