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).
I have written about this model previously in this stan help forum here: Where to declare individual survival parameter in multi-state state-space model (Again, many thanks to Daniel Lee for his quick reply.)
My model is working as expected when the response variable is completely observed.
However, I have the perhaps slightly unusual case that observations at one of the occasions are missing because a detection device was unavailable.
I have read chapter 10 (page 177+ in stan reference manual 2.15.0 ) and I have searched this forum and studied the following hits:
But still I have been unable to find a solution to my issue.
Here is my (simplified) model:
// -------------------------------------------------
// 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];
}
// Transformed data
transformed data {
// Varaible calculations
int n_occs_less1 = n_occs - 1;
}
// Parameters
parameters {
// Parameter declarations
real<lower = 0, upper = 1> phi_juv; // ad survival
real<lower = 0, upper = 1> phi_ad; // ad survival
real<lower = 0, upper = 1> p; // probability of detection
}
// Transformed parameters
transformed parameters {
// Parameter declarations
simplex[3] ps[3, n_inds, n_occs_less1];
simplex[3] po[3, n_inds, n_occs_less1];
// 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;
ps[1, i, o, 3] = (1 - phi_juv);
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];
// 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
And here is my model with attempts to account for the missing data at occasion 3 in the first year:
// -------------------------------------------------
// 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 = 0> n_obscol3;
int<lower = 0> n_miscol3;
int<lower = 1, upper = n_obscol3 + n_miscol3> i_obscol3[n_obscol3];
int<lower = 1, upper = n_obscol3 + n_miscol3> i_miscol3[n_miscol3];
vector<lower = 1, upper = 3>[n_obscol3] y_obscol3;
int<lower = 1, upper = 3> y_obs[n_inds, 2];
}
// Transformed data
transformed data {
// Varaible calculations
int n_occs_less1 = n_occs - 1;
}
// Parameters
parameters {
// Parameter declarations
real<lower = 0, upper = 1> phi_juv; // juv survival
real<lower = 0, upper = 1> phi_ad; // ad survival
real<lower = 0, upper = 1> p; // probability of detection
vector<lower = 1, upper = 3>[n_miscol3] y_miscol3;
}
// Transformed parameters
transformed parameters {
// Parameter declarations
simplex[3] ps[3, n_inds, n_occs_less1];
simplex[3] po[3, n_inds, n_occs_less1];
// matrix<lower = 1, upper = 3>[n_inds, n_occs] y;
int<lower = 1, upper = 3> y[n_inds, n_occs];
y[i_obscol3, 3] = y_obscol3;
y[i_miscol3, 3] = y_miscol3;
// 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;
ps[1, i, o, 3] = (1 - phi_juv);
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];
// 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
I am getting the following error during compilation:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
parameters or transformed parameters cannot be integer or integer array; found declared type int, parameter name=y
Problem with declaration.
error in 'model35d018f9625_eebc87a8ea17142b0af29d6ae65307d5' at line 55, column 48
-------------------------------------------------
53: simplex[3] po[3, n_inds, n_occs_less1];
54: // matrix<lower = 1, upper = 3>[n_inds, n_occs] y;
55: int<lower = 1, upper = 3> y[n_inds, n_occs];
^
56: y[i_obscol3, 3] = y_obscol3;
-------------------------------------------------
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model 'eebc87a8ea17142b0af29d6ae65307d5' due to the above error.
To me, this suggests that the following definition is not allowed in the transformed parameters or parameters block:
int<lower = 1, upper = 3> y[n_inds, n_occs];
When I replace it with:
matrix<lower = 1, upper = 3>[n_inds, n_occs] y;
in the transformed parameter block, the compilation proceeds into the model block but then I get the error that y should be an integer:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
ERROR: Container index must be integer; found type=real
error in 'model35d06ed8176b_42b1f38c990210b6cf618b2cb4d895a2' at line 112, column 77
-------------------------------------------------
110: for (k in 1:3) {
111: for (j in 1:3) {
112: acc[j] = gamma[o - 1, j] * ps[j, i, o - 1, k] * po[k, i, o - 1, y[i, o]];
^
113: }
-------------------------------------------------
PARSER EXPECTED: <one or more container indexes followed by ']'>
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model '42b1f38c990210b6cf618b2cb4d895a2' due to the above error.
which makes sense.
Here is some R code to simulate some data:
## define function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, marked, nocc){
nyr <- length(marked)
CH <- CH.TRUE <- array(NA, dim = list(max(marked), nocc, nyr))
for(y in 1:nyr){
CH[1:marked[y], 1, y] <- 1 # always seen at release
CH.TRUE[1:marked[y], 1, y] <- 1 # always seen at release
}
for(y in 1:nyr){
for(i in 1:marked[y]){
for(o in 2:nocc){
state.p <- rmultinom(1, 1, PSI.STATE[[y]][CH.TRUE[i, o-1, y],,i,o-1])
state <- which(state.p == 1)
CH.TRUE[i, o, y] <- state
event.p <- rmultinom(1, 1, PSI.OBS[[y]][CH.TRUE[i, o, y],,i,o-1])
event <- which(event.p == 1)
CH[i, o, y] <- event
} #o
} #i
} #y
return(list(CH = CH, CH.TRUE = CH.TRUE))
}
## define number of states, occasions, released individuals & years
n.states <- 3 # juv, ad, dead
n.occs <- 3
n.years <- 3
marked <- sample(x = 100:300, size = n.years, replace = TRUE)
years <- factor(rep(1:n.years, marked), labels = (2017 - (n.years - 1)):2017)
n.inds <- sum(marked)
## define parameters
phi.juv <- 0.08
phi.ad <- 0.8
p <- 0.75
## define state and observation matrices
PSI.STATE <- vector('list', n.years)
for(y in 1:n.years){
PSI.STATE[[y]] <- array(NA, dim = c(n.states, n.states, marked[y], n.occs - 1))
}
for(y in 1:n.years){
for(i in 1:marked[y]){
for(o in 1:(n.occs - 1)){
PSI.STATE[[y]][,,i,o] <- matrix(c(
0, phi.juv, 1 - phi.juv,
0, phi.ad, 1 - phi.ad,
0, 0, 1
), nrow = n.states, byrow = TRUE)
} # o
} # i
} # y
PSI.OBS <- vector('list', n.years)
for(y in 1:n.years){
PSI.OBS[[y]] <- array(NA, dim = c(n.states, n.states, marked[y], n.occs - 1))
}
for(y in 1:n.years){
for(i in 1:marked[y]){
for(o in 1:(n.occs - 1)){
PSI.OBS[[y]][,,i,o] <- matrix(c(
0, 0, 1,
0, p, 1-p,
0, 0, 1
), nrow = n.states, byrow = TRUE)
} # o
} # i
} # y
## execute simulation function
sim <- simul.ms(PSI.STATE, PSI.OBS, marked, n.occs)
ch <- sim$CH
## make into a single matrix
ch.mat <- matrix(NA, ncol = (n.years + (n.occs - 1)), nrow = sum(marked))
ridx <- 1
cidx <- 1
for(y in 1:n.years){
ch.mat[ridx:(ridx + (marked[y] - 1)), cidx:(cidx + (n.occs - 1))] <- ch[1:marked[y], , cidx]
ridx <- ridx + marked[y]
cidx <- cidx + 1
}
## recode NAs as state 3
ch.mat[is.na(ch.mat)] <- 3
## add NAs at occasion 3 in year 1
if(addNAs){
ch[, n.occs, 1] <- NA
}
## make an easy matrix
ch.emat <- numeric()
for(i in 1:dim(ch)[3]){
ch.emat <- rbind(ch.emat, ch[1:marked[i], , i])
}
## dice it up into constituent parts
y.obs <- ch.emat[, 1:2]
i.obscol3 <- which(!is.na(ch.emat[, 3]))
i.miscol3 <- which(!is.na(ch.emat[, 3]))
n.obscol3 <- length(i.obscol3)
n.miscol3 <- length(i.miscol3)
y.obscol3 <- as.integer(na.omit(ch.emat[, 3]))
Could anyone help / guide me towards a possible solution to this problem?
Many thanks in advance for your help,
Stephen
p.s., sorry for the obscenely long post!