Where to declare individual survival parameter in multi-state state-space model

Hello,

I am trying to fit a multi-state mark-recapture state-space model, similar to those in Chapter 9 of Marc Kery’s BPA (translated to stan code at https://github.com/stan-dev/example-models/tree/master/BPA) but including an individual continuous covariate (length).

My model is as follows

  // -------------------------------------------------
  // States (S):
  // 1 juvenile
  // 2 adult
  // 3 dead
  // Observations (O):
  // 1 seen as juvenile
  // 2 seen as adult
  // 3 not seen
  // -------------------------------------------------

  // Data
  data {
    
    // Variable declarations
    int<lower = 0> n_yrs;
    int<lower = 0> n_inds;
    int<lower = 0> n_occs;
    int<lower = 1, upper = 3> y[n_inds, n_occs];
    int years[n_inds];
    real lengths[n_inds];

  }

  // Transformed data
  transformed data {
  
    // Varaible calculations
    int n_occs_less1 = n_occs - 1;
    
  }

  // Parameters
  parameters {
    
    // Parameter declarations
    real<lower = 0, upper = 1> phi_ad;             // ad survival
    real<lower = 0, upper = 1> p;                  // probability of detection
    real a1[n_yrs];
    real a1_mu;
    real<lower = 0> a1_sigma;
    real b1;
    // real<lower = 0, upper = 1> phi_juv[n_inds];    // juv survival

  }

  // Transformed parameters
  transformed parameters {

    // Parameter declarations
    simplex[3] ps[3, n_inds, n_occs_less1];
    simplex[3] po[3, n_inds, n_occs_less1];
    real<lower = 0, upper = 1> phi_juv[n_inds];    // juv survival
    real lp[n_inds];
    
    // Constraints
    for (i in 1:n_inds) {
      phi_juv[i] = inv_logit(lp[i]);
      lp[i] = a1[years[i]] + (b1 * lengths[i]);
    }
    print(phi_juv[1:n_inds]);
    print(lp[1:n_inds]);

    // Define state-transition and observation matrices
    for (i in 1:n_inds) {
      for (o in 1:n_occs_less1) {
        
        // Define probabilities of state S(y+1) given S(y)
        ps[1, i, o, 1] = 0.0;
        // ps[1, i, o, 2] = phi_juv[i];
        // ps[1, i, o, 3] = (1 - phi_juv[i]);
        ps[1, i, o, 2] = 0.5;
        ps[1, i, o, 3] = 0.5;
        ps[2, i, o, 1] = 0.0;
        ps[2, i, o, 2] = phi_ad;
        ps[2, i, o, 3] = (1 - phi_ad);
        ps[3, i, o, 1] = 0.0;
        ps[3, i, o, 2] = 0.0;
        ps[3, i, o, 3] = 1.0;

        // Define probabilities of O(y) given S(y)
        po[1, i, o, 1] = 0.0;
        po[1, i, o, 2] = 0.0;
        po[1, i, o, 3] = 1.0;
        po[2, i, o, 1] = 0.0;
        po[2, i, o, 2] = p;
        po[2, i, o, 3] = (1 - p);
        po[3, i, o, 1] = 0.0;
        po[3, i, o, 2] = 0.0;
        po[3, i, o, 3] = 1.0;

      } // o
    } // i

  }

  // Model
  model {

    // Variable declarations
    real acc[3];
    vector[3] gamma[n_occs];

    // Priors
    for (t in 1:n_yrs) {
      a1[t] ~ normal(a1_mu, a1_sigma);
    }
    a1_mu ~ normal(0, 100);
    a1_sigma ~ uniform(0, 100);
    b1 ~ normal(0, 10);
    phi_ad ~ uniform(0, 1);
    p ~ uniform(0, 1);

    // Likelihood
    // Forward algorithm derived from Stan Modeling Language User's Guide and Reference Manual
    for (i in 1:n_inds) {
      for (k in 1:3) {
        gamma[1, k] = (k == y[i, 1]);
      }
      for (o in 2:n_occs) {
        for (k in 1:3) {
          for (j in 1:3) {
            acc[j] = gamma[o - 1, j] * ps[j, i, o - 1, k] * po[k, i, o - 1, y[i, o]];
          }
          gamma[o, k] = sum(acc);
        }
      }
      target += log(sum(gamma[n_occs]));
    }

  } // end model

The code compiles but the model fails to initialise (with error: “Exception: validate transformed params: phi_juv[k0__] is nan, but must be greater than or equal to 0”), I think because phi_juv is declared in the Transformed Parameters block rather than the Parameters block. However, I can’t see how I can declare it in the Parameter block and still apply the constraint that it is a logit of covariates.

Can anyone help?
I would really appreciate some guidance as I am new to Stan but would like to use it.

Many thanks in advance for your help,
Stephen

The Stan language is imperative e.g. every line is executed in order.

I think you’ll get around your first problem by swapping these two lines:

Want to try that first?

Hi syclik,

Thanks for your super quick reply.
Sorry for wasting your time(!) - it seems to work perfectly now.

Perhaps you can tell that I am coming from a BUGS background.

Many thanks again,
Stephen

1 Like

Mind marking my reply as the solution?

You’re welcome! It wasn’t a waste of time.

There were a few things in your post that made it easy to answer on the first read:

  1. Described what you’re trying to do at a high-level
  2. Included the full Stan program in a code block (when it’s embedded in a file, it’s yet another click and another barrier to answering a question quickly)
  3. Included the full error message; it’s a lot easier for people to help when the error message is included!
  4. A guess as to what you thought could be wrong

Thanks for making it super-easy to answer.

=).

1 Like

Still very quick!

I marked your answer as the solution.

Many thanks again!

1 Like

This is all really key. We say it in web pages on how to report issues and so on, but it’s hard to get the instructions in front of people at the right time!

1 Like