Hi, I am attempting to roughly translate a JAGS model from Olivier Gimenez’s Github page, Bayesian implementation of capture-recapture models with robust design. Here is the JAGS code:
model <- function() {
# priors
phi ~ dunif(0,1) # survival
gP ~ dunif(0,1) # gamma'
gPP ~ dunif(0,1) # gamma''
gamP <- 1 - gP # MARK parameterization
gamPP <- 1 - gPP # MARK parameterization
mean.p ~ dunif(0,1) # detection
# secondary occasions p's
for (t in 1:n.years){
for (j in 1:max(n.sec[1:n.years])){
p[t,j] <- mean.p
}
}
for (t in 1:n.years){
for (j in 1:n.sec[t]){
yes[t,j] ~ dbin(p[t,j], total[t,j])
}
}
# Primary occasions p's or pooled detection probability
for (t in 1:n.years){
pstar[t] <- 1 - prod(1 - p[t,1:n.sec[t]])
}
# state matrices
s[1,1] <- phi * gPP
s[1,2] <- phi * (1 - gPP)
s[1,3] <- 1 - phi
s[2,1] <- phi * gP
s[2,2] <- phi * (1 - gP)
s[2,3] <- 1 - phi
s[3,1] <- 0
s[3,2] <- 0
s[3,3] <- 1
# observation matrices
for (t in 1:n.years){
o[1,t,1] <- pstar[t]
o[1,t,2] <- 1 - pstar[t]
o[2,t,1] <- 0
o[2,t,2] <- 1
o[3,t,1] <- 0
o[3,t,2] <- 1
}
# likelihood
for (i in 1:n.ind){
z[i,first[i]] <- ch[i,first[i]]
for (t in (first[i]+1):n.years){
z[i,t] ~ dcat(s[z[i,t-1], ]) # state equations
ch[i,t] ~ dcat(o[z[i,t], t, ]) # obsevation equations
}
}
}
My main issue is that I receive the error, “…evaluating the log probability at the initial value…” because ‘p’ is nan. Somehow I need to calculate p for secondary capture occasions first and then evaluate the primary capture occasions in the rest of the model but I do not know where to go with the binomial
for (t in 1:n.years){
for (j in 1:n.sec[t]){
yes[t,j] ~ dbin(p[t,j], total[t,j])
}
}
part of the JAGS model so that the transformed parameters then have the recapture probability ‘p’ to work with.
// ---------------------------------
// States (S):
// 1 alive and present
// 2 alive and absent
// 3 dead
// Observations (O):
// 1 seen
// 2 not seen
// ---------------------------------
functions {
/**
* Return an integer value denoting occasion of first capture.
* This function is derived from Stan Modeling Language
* User's Guide and Reference Manual.
*
* @param y Observed values
* @return Occasion of first capture
*/
int first_capture(array[] int y_i) {
for (k in 1 : size(y_i)) {
if (y_i[k] == 1) {
return k;
}
}
return 0;
}
}
data {
int<lower=0> nind; //number of individuals
//int<lower=0> n_occasions; //number of occassions
int<lower=1> nprimary; // Number of years or primary occasions in this case
vector<lower=1>[nprimary] nsecondary; //number of secondary occassions
int<lower=1> nsecondary_max; //max number of secondary occasions
array[nind, nprimary] int<lower=1,upper=2> ch; //individual capture histories for each primary occasion
array[nprimary, nsecondary_max] int<lower=0> yes; //fish recaptured by primary and secondary occasion
array[nprimary, nsecondary_max] int<lower=0> total; //total fish caught during primary and secondary occasion
}
transformed data {
//int n_occ_minus_1 = n_occasions - 1;
int<lower=0,upper=nprimary> first[nind];
for (i in 1:nind)
first[i] = first_capture(ch[i]);
}
parameters {
real<lower=0,upper=1> mean_phi; // mean survival
real<lower=0,upper=1> mean_gP; // mean gamma'
real<lower=0,upper=1> mean_gPP; // mean gamma''
real<lower=0,upper=1> mean_p; // mean recapture probability
}
transformed parameters {
vector<lower=0,upper=1>[nprimary] phi; //survival probability
vector<lower=0,upper=1>[nprimary] gP; // the probability of being off the study area, unavailable for capture during primary trapping session (i) given that the animal was not present on the study area during primary trapping session (8 ??? 1), and survives to trapping session (8).
vector<lower=0,upper=1>[nprimary] gPP; // the probability of being off the study area, unavailable for capture during the primary trapping session (i) given that the animal was present during primary trapping session (8 ??? 1), and survives to trapping session (8).
matrix<lower=0,upper=1>[nind, nsecondary_max] p; // Recapture probability
vector<lower=0,upper=1>[nprimary] pstar; //primary occasion p's or pooled detection probability
array[3, nprimary] simplex[3] ps;
array[3, nprimary] simplex[2] po;
// Constraints
for (t in 1 : nprimary) {
phi[t] = mean_phi;
gP[t] = mean_gP;
gPP[t] = mean_gPP;
}
//secondary occasions p's
for (t in 1:nprimary){
for (j in 1:nsecondary_max){
p[t,j] = mean_p;
}
}
//Primary occasions p's or pooled detection probability
for (t in 1:nprimary){
pstar[t] = 1 - prod(1 - p[t,1:nsecondary_max]);
}
// Define state-transition and observation matrices
for (t in 1 : nprimary) {
// Define probabilities of state S(t+1) given S(t)
ps[1, t, 1] = phi[t] * gPP[t];
ps[1, t, 2] = phi[t] * (1 - gPP[t]);
ps[1, t, 3] = 1.0 - phi[t];
ps[2, t, 1] = phi[t] * gP[t];
ps[2, t, 2] = phi[t] * (1 - gP[t]);
ps[2, t, 3] = 1.0 - phi[t];
ps[3, t, 1] = 0;
ps[3, t, 2] = 0;
ps[3, t, 3] = 1.0;
// Define probabilities of O(t) given S(t)
po[1, t, 1] = pstar[t];
po[1, t, 2] = 1 - pstar[t];
po[2, t, 1] = 0.0;
po[2, t, 2] = 1.0;
po[3, t, 1] = 0.0;
po[3, t, 2] = 1.0;
}
}
model {
array[nind, nprimary] int z;
// priors
// phi ~ uniform(0,1); # survival
// gP ~ uniform(0,1); # gamma'
// gPP ~ uniform(0,1); # gamma''
// mean_p ~ uniform(0,1); # detection
//Likelihood for secondary occasion p's
for (t in 1:nprimary){
for (j in 1:nsecondary_max){
yes[t,j] ~ binomial(total[t,j], p[t,j]);
}
}
// Likelihood
for (i in 1:nind){
z[i,first[i]] = ch[i,first[i]];
for (t in (first[i]+1):nprimary){
z[i,t] ~ categorical(ps[i,t-1]); # state equations
ch[i,t] ~ categorical(po[i,t]); # obsevation equations
}
}
}