insufficient description of input data for individual CJS mark-recapture model in stan user guide

47 Views Asked by At

EDIT This issue appears to have resulted from rstan not having up-to-date stan syntax because of the restrictions that CRAN places on package updates. A more up-to-date rstan package appears to be available at the mc-stan repo, but does not yet have binaries for the most recent R version (4.3).


I am not sure exactly where this question goes, but haven't found an answer anywhere else, so here goes.

I need to learn stan to run some CJS mark-recapture models, and have been working from the stan user guide here to do so. The collective Cormack-Jolly-Seber model works fine, and is fairly straightforward in terms of the data going into it, but I really need to get a version of an individual CJS working. The problem is that the guide does a very poor job of describing the data that goes into the model and there is no data provided. As a result, when I run the example model with real data, it keeps throwing back errors that suggest something is wrong with the input data.

Here is the data block, which requires a capture history matrix y of I rows and T columns:

data {
  int<lower=2> T;
  int<lower=0> I;
  int<lower=0,upper=1> y[I, T];
}

From the parameters, I know that it also requires phi and p vectors of lengths T-1 and T, respectively.

parameters {
  vector<lower=0,upper=1>[T-1] phi;
  vector<lower=0,upper=1>[T] p;
}

So far this is consistent with the collective CJS model and shouldn't be causing issues. I have also included values for T and I because of their inclusion in the data block. So so far my input data looks like this (I am using rstan, so this is a toy version of what it looks like in r):

history.array<-list(y=matrix(c(1,0,0,1,0,0,0,0,1,0,0,1,
         0,0,0,0,1,1,0,1,1,0,0,1,
         1,1,0,0,1,1,0,0,0,0,0,0,
         1,1,0,0,1,0,0,1), nrow=11),
phi=rep(.5,4-1),
p=rep(.5,4),
T = as.integer(4),
I = as.integer(11))

And here is the full of the stan code from the guide (which I have saved as "individual_CJS.stan").

data {
  int<lower=2> T;
  int<lower=0> I;
  int<lower=0, upper=1> y[I,T];
}

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;
      k = size(y_i) - k_rev;
      if (y_i[k])
        return k;
    }
    return 0;
  }
  vector prob_uncaptured(int T, vector p, vector phi) {
    vector[T] chi;
    chi[T] = 1.0;
    for (t in 1:(T - 1)) {
      int t_curr;
      int t_next;
      t_curr = T - t;
      t_next = t_curr + 1;
      chi[t_curr] = (1 - phi[t_curr])
                     + phi[t_curr]
                       * (1 - p[t_next])
                       * chi[t_next];
    }
    return chi;
  }
}

transformed data {
  int<lower=0,upper=T> first[I];
  int<lower=0,upper=T> last[I];
  vector<lower=0,upper=I>[T] n_captured;
  for (i in 1:I)
    first[i] = first_capture(y[i]);
  for (i in 1:I)
    last[i] = last_capture(y[i]);
  n_captured = rep_vector(0, T);
  for (t in 1:T)
    for (i in 1:I)
      if (y[i, t])
        n_captured[t] = n_captured[t] + 1;
}

parameters {
  vector<lower=0,upper=1>[T-1] phi;
  vector<lower=0,upper=1>[T] p;
}
transformed parameters {
  vector<lower=0,upper=1>[T] chi;
  chi = prob_uncaptured(T,p,phi);
}

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

generated quantities {
  real beta;
  vector<lower=0>[T] pop;

  beta = phi[T - 1] * p[T];
  pop = n_captured ./ p;
  pop[1] = -1;
}

But when I try to run the program with these inputs:

warmups <- 1000
total_iterations <- 2000
max_treedepth <-  10
n_chains <-  4
n_cores <- 4
adapt_delta <- 0.95
test_ind_CJS<-stan(
  file = "individual_CJS.stan",
  chains = n_chains,
  warmup = warmups,
  data = history.array,
  iter = total_iterations,
  cores = n_cores,
  refresh = 250,
  control = list(max_treedepth = max_treedepth,
                 adapt_delta = adapt_delta)
),

I get this error:

PARSER FAILED TO PARSE INPUT COMPLETELY
STOPPED AT LINE 10: 
functions {
  int first_capture(int[] y_i) {
    for (k in 1:size(y_i))
      if (y_i[k])
        return k;
    return 0;
  }.......*the rest of the stan code is printed here*.......

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'individual_CJS' due to the above error.

So that first_capture function seems to need something, but I can't figure out what. y_i, shouldn't need to be input if it is just the ith row of matrix y. What am I missing here? thanks in advance

0

There are 0 best solutions below