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