Error when running multistate Jolly-Seber model from BPA

[Matthijs Hollanders]

Hi all,

I'm trying to run a multistate Jolly-Seber model in JAGS using Kéry & Schaub BPA. To start, I simply copied the BPA JAGS code, specifically 1) the model under section 10.3.2, 2) the data generation code under section 10.4, and 3) the code to run the model under 10.4.2. 

When I run this code, I get the following error after running JAGS:

"Error in setParameters(init.values[[i]], i) : Error in node po[1,1,1,2]

Cannot set value of non-variable node"

FYI, when I do the same for the restricted occupancy and superpopulation parameterisations, the models run just fine. I do actually need to use the multistate formulation because the model will later be expanded to include multiple disease states. Can anyone clarify what's going on? 

Looking forward to your answers. :)

Kind regards,

Matt Hollanders

[Matthijs Hollanders]

Since it's not letting me attach the script, here's the code I've used copied:

library(jagsUI)

# 10.3.2. The JS model as a multistate model

# Specify model in BUGS language

sink("js-ms.jags")

cat("

model {

#--------------------------------------

# Parameters:

# phi: survival probability

# gamma: removal entry probability

# p: capture probability

#--------------------------------------

# States (S):

# 1 not yet entered

# 2 alive

# 3 dead

# Observations (O):

# 1 seen 

# 2 not seen

#--------------------------------------

# Priors and constraints

for (t in 1:(n.occasions-1)){

   phi[t] <- mean.phi

   gamma[t] ~ dunif(0, 1) # Prior for entry probabilities

   p[t] <- mean.p

   }

mean.phi ~ dunif(0, 1)    # Prior for mean survival

mean.p ~ dunif(0, 1)      # Prior for mean capture

# Define state-transition and observation matrices

for (i in 1:M){  

   # Define probabilities of state S(t+1) given S(t)

   for (t in 1:(n.occasions-1)){

      ps[1,i,t,1] <- 1-gamma[t]

      ps[1,i,t,2] <- gamma[t]

      ps[1,i,t,3] <- 0

      ps[2,i,t,1] <- 0

      ps[2,i,t,2] <- phi[t]

      ps[2,i,t,3] <- 1-phi[t]

      ps[3,i,t,1] <- 0

      ps[3,i,t,2] <- 0

      ps[3,i,t,3] <- 1

      

      # Define probabilities of O(t) given S(t)

      po[1,i,t,1] <- 0

      po[1,i,t,2] <- 1

      po[2,i,t,1] <- p[t]

      po[2,i,t,2] <- 1-p[t]

      po[3,i,t,1] <- 0

      po[3,i,t,2] <- 1

      } #t

   } #i

# Likelihood 

for (i in 1:M){

   # Define latent state at first occasion

   z[i,1] <- 1   # Make sure that all M individuals are in state 1 at t=1

   for (t in 2:n.occasions){

      # State process: draw S(t) given S(t-1)

      z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])

      # Observation process: draw O(t) given S(t)

      y[i,t] ~ dcat(po[z[i,t], i, t-1,])

      } #t

   } #i

# Calculate derived population parameters

for (t in 1:(n.occasions-1)){

   qgamma[t] <- 1-gamma[t]

   }

cprob[1] <- gamma[1]

for (t in 2:(n.occasions-1)){

   cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)])

   } #t

psi <- sum(cprob[])            # Inclusion probability

for (t in 1:(n.occasions-1)){

   b[t] <- cprob[t] / psi      # Entry probability

   } #t

for (i in 1:M){

   for (t in 2:n.occasions){

      al[i,t-1] <- equals(z[i,t], 2)

      } #t

   for (t in 1:(n.occasions-1)){

      d[i,t] <- equals(z[i,t]-al[i,t],0)

      } #t   

   alive[i] <- sum(al[i,])

   } #i

for (t in 1:(n.occasions-1)){

   N[t] <- sum(al[,t])        # Actual population size

   B[t] <- sum(d[,t])         # Number of entries

   } #t

for (i in 1:M){

   w[i] <- 1-equals(alive[i],0)

   } #i

Nsuper <- sum(w[])            # Superpopulation size

}

",fill = TRUE)

sink()

# 10.4. Models with constant survival and time-dependent entry

# Define parameter values

n.occasions <- 7                         # Number of capture occasions

N <- 400                                 # Superpopulation size

phi <- rep(0.7, n.occasions-1)           # Survival probabilities

b <- c(0.34, rep(0.11, n.occasions-1))   # Entry probabilities 

p <- rep(0.5, n.occasions)               # Capture probabilities

PHI <- matrix(rep(phi, (n.occasions-1)*N), ncol = n.occasions-1, nrow = N, byrow = T)

P <- matrix(rep(p, n.occasions*N), ncol = n.occasions, nrow = N, byrow = T)

# Function to simulate capture-recapture data under the JS model

simul.js <- function(PHI, P, b, N){

  B <- rmultinom(1, N, b) # Generate no. of entering ind. per occasion

  n.occasions <- dim(PHI)[2] + 1

  CH.sur <- CH.p <- matrix(0, ncol = n.occasions, nrow = N)

  # Define a vector with the occasion of entering the population

  ent.occ <- numeric()

  for (t in 1:n.occasions){

    ent.occ <- c(ent.occ, rep(t, B[t]))

  }

  # Simulating survival

  for (i in 1:N){

    CH.sur[i, ent.occ[i]] <- 1   # Write 1 when ind. enters the pop.

    if (ent.occ[i] == n.occasions) next

    for (t in (ent.occ[i]+1):n.occasions){

      # Bernoulli trial: has individual survived occasion?

      sur <- rbinom(1, 1, PHI[i,t-1])

      ifelse (sur==1, CH.sur[i,t] <- 1, break)

    } #t

  } #i

  # Simulating capture

  for (i in 1:N){

    CH.p[i,] <- rbinom(n.occasions, 1, P[i,])

  } #i

  # Full capture-recapture matrix

  CH <- CH.sur * CH.p

  

  # Remove individuals never captured

  cap.sum <- rowSums(CH)

  never <- which(cap.sum == 0)

  CH <- CH[-never,]

  Nt <- colSums(CH.sur)    # Actual population size

  return(list(CH=CH, B=B, N=Nt))

}

# Execute simulation function

sim <- simul.js(PHI, P, b, N)

CH <- sim$CH

# 10.4.2 Analysis of the JS model as a multistate model

# Add dummy occasion

CH.du <- cbind(rep(0, dim(CH)[1]), CH)

# Augment data

nz <- 500

CH.ms <- rbind(CH.du, matrix(0, ncol = dim(CH.du)[2], nrow = nz))

# Recode CH matrix: a 0 is not allowed in WinBUGS!

CH.ms[CH.ms==0] <- 2                     # Not seen = 2, seen = 1

# Bundle data

jags.data <- list(y = CH.ms, n.occasions = dim(CH.ms)[2], M = dim(CH.ms)[1])

# Initial values

# As always in JAGS, good initial values need to be specified for the latent state z. They need to correspond to the true state, which is not the same as the observed state. Thus, we have given initial values of "1" for the latend state at all places before an individual was observed, an initial value of "2" at all places when the individual was observed alive or known to be alive and an initial value of "3" at all places after the last observation. The folllowing function creates the initial values.

js.multistate.init <- function(ch, nz){

  ch[ch==2] <- NA

  state <- ch

  for (i in 1:nrow(ch)){

    n1 <- min(which(ch[i,]==1))

    n2 <- max(which(ch[i,]==1))

    state[i,n1:n2] <- 2

  }

  state[state==0] <- NA

  get.first <- function(x) min(which(!is.na(x)))

  get.last <- function(x) max(which(!is.na(x)))   

  f <- apply(state, 1, get.first)

  l <- apply(state, 1, get.last)

  for (i in 1:nrow(ch)){

    state[i,1:f[i]] <- 1

    if(l[i]!=ncol(ch)) state[i, (l[i]+1):ncol(ch)] <- 3

    state[i, f[i]] <- 2

  }   

  state <- rbind(state, matrix(1, ncol = ncol(ch), nrow = nz))

  return(state)

}

inits <- function(){list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = js.multistate.init(CH.du, nz))}    

# Parameters monitored

parameters <- c("mean.p", "mean.phi", "b", "Nsuper", "N", "B")

# MCMC settings

ni <- 20000

nt <- 3

nb <- 5000

nc <- 3

# Call JAGS from R (BRT 32 min)

js.ms <- jags(jags.data, inits, parameters, "js-ms.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)

[macedo...@gmail.com]

Hi Matt,

I get the same error but I'm trying to run a Dail-Madsen model. Did you find out how to solve this problem?

Kind regards,

Thais

[Jose Jimenez Garcia-Herrera]

Hi Matt,

 

I don’t know is too later for you, but here you are the code working correctly. Maybe it can be useful for others.

You can find the original solution for a similar issue at this page:

 

https://stackoverflow.com/questions/45202871/jolly-seber-runs-with-winbugs-not-with-jags

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2020)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hmecology/b55e6471-b09c-4801-94f0-37bc0c73d77an%40googlegroups.com.