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
// Observations (O):
// 1 seen as juvenile
// 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> 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);
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.

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,

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.

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!