Uneven time intervals in mark-recapture

Modeling
,
1

Dear All,
I am new to Stan and I am using the RStan models from Kéry and Schaub 2012 as kindly provided by Hiroki Itô. I am currently using a time-varying CJS model, and would like to adjust my quarterly (3 month) survival estimates to account for an uneven sampling event of 2 months followed by a period of 4 months.
I think I want to estimate a monthly survival, which can then be multiplied to calculate survival for the 2 month and 4 month periods.
phi_month = exp(log(3_month_phi)/3)
phi_2_months = phi_month^2
phi_4_months = phi_month^4

There are also a few incidences of missed sampling.

The unaltered model.

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)) {
      // Compoud declaration was enabled in Stan 2.13
      int k = size(y_i) - k_rev;
      //      int k;
      //      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;
        /*
        int t_curr;
        int t_next;

        t_curr = n_occasions - t;
        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];
  int<lower=0,upper=1> p_deploy[24];
  }

transformed data {
  int n_occ_minus_1;
  // Compoud declaration is enabled in Stan 2.13
  //  int n_occ_minus_1 = n_occasions - 1;
  int<lower=0,upper=n_occasions> first[nind];
  int<lower=0,upper=n_occasions> last[nind];

  n_occ_minus_1 = n_occasions - 1;
  for (i in 1:nind)
    first[i] = first_capture(y[i]);
  for (i in 1:nind)
    last[i] = last_capture(y[i]);
}

parameters {
  vector<lower=0,upper=1> [n_occ_minus_1] alpha;   // Mean survival
  vector<lower=0,upper=1>[n_occ_minus_1] beta;   // Mean recapture
}

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;

  // 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] = alpha[t];
      p[i, t] = beta[t];
    }
  }

  chi = prob_uncaptured(nind, n_occasions, p, phi);
}

model {
  // Priors
  // Uniform priors are implicitly defined.
  //  alpha ~ uniform(0, 1);
  //  beta ~ uniform(0, 1);

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

  }    
}

I’m uncertain of my solution which is to include a line in the constraints, something like…

  // 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_int[t]*(log(alpha[t])/3)); // phi_int[i, t] is a vector of time intervals e.g. 3,3,3,2,4,3,,,
      p[i, t] = beta[t]*p_deploy[t]; // p_deploy is a vector of 1s and 0s as there is some missed sampling
    }
  }

Any help forum members could offer that would help me clarify my thinking on this would be very much appreciated.

2

Since nobody else answered, I will give it a try:

I am not at all familiar with those models (or actually probably your entire field, sorry :-) ) providing some background/links/… on those models is a good idea, as the crowd here tends to be very diverse and few people share the same background.

I however think I was able to understand what the model does from the code. My best guess (though I am no expert) is that your adjustments seem reasonable. The best way to check is (as always) to write code to simulate data where you know the alpha and beta values and see if the model recovers them correctly.

Also, I am quite sure that you can simplify phi[i, t] = exp(phi_int[t]*(log(alpha[t])/3)); to phi[i, t] = alpha_t ^ (to_real(phi_int) / 3); (but please double check my math).

Hope that helps.

3

Thank you for the reassurance, and for the suggested simplification. It’s really helpful to get feedback, I’ll try and provide a bit more context next time. Thanks again!