Creating a derived parameter from sims.list within jags model

Question

I am working on a combined integrated model/population viability analysis in jags, and am trying to create a derived parameter to use in my population projections.

In the current version of the model that runs without errors, I have time-dependent survival parameters that are generated within the capture-mark-recapture portion of the model, and feed into the population state-space model for the years I have data for (t in 1:n.occasions-1). I have created derived parameters based on the mean of these parameters, that I have used in my population state-space model for the years I would like to project over (t in n.occasion:n.occasions-1+K). However, by using this approach, I think the variance around the mean of the derived parameters is too small and is affecting my predicted population sizes.

As such, I would like to try creating a derived parameter that is a random draw from the sims.list for the time-dependent survival parameters, to capture the variance better. My question essentially, is can you draw from the sims list of a model, within the model itself? Or can you only draw from that sims list once the model has finished running?

In terms of what I have tried and the errors I receive:

I have tried creating a derived parameter, and get a syntax error near "$", relating to the line where I reference the sims.list.

My model is over 3000 lines long, so I have included an example of a simpler model below, where I have tried to do the same thing, and receive the same error (Error in checkForRemoteErrors(val) : 3 nodes produced errors; first error: Error parsing model file: syntax error on line 103 near "$), to illustrate the issue I am facing. The model code I am using is courtesy of Michael Schaub and Marc Kery from their IPM Book. In this example, I have created the derived parameter sj.f, which I am using to replace sj[t] in the part of the population model that predicts into the future.

This is my first question on stackoverflow (and second attempt at posting, this time with a reproducible example), so apologies if I have included too little or too much information.

# Load woodchat shrike data and produce data overview
library(IPMbook)
data(woodchat10)
str(woodchat10)

K <- 15    # Number of forecast years
jags.data <- list(marr.j=woodchat10$marr.j, marr.a=woodchat10$marr.a, n.occasions=ncol(woodchat10$marr.j), rel.j=rowSums(woodchat10$marr.j), rel.a=rowSums(woodchat10$marr.a), J=woodchat10$J, B=woodchat10$B, count=woodchat10$count, pNinit=dUnif(1, 50), K=K)
str(jags.data)


K <- 15    # Number of forecast years
jags.data <- list(marr.j=woodchat10$marr.j, marr.a=woodchat10$marr.a, n.occasions=ncol(woodchat10$marr.j), rel.j=rowSums(woodchat10$marr.j), rel.a=rowSums(woodchat10$marr.a), J=woodchat10$J, B=woodchat10$B, count=woodchat10$count, pNinit=dUnif(1, 50), K=K)
str(jags.data)


# Write JAGS model file
cat(file="model1.txt", "
model { 
# Priors and linear models
mean.logit.sj <- logit(mean.sj)
mean.sj ~ dunif(0, 1) 
mean.logit.sa <- logit(mean.sa)
mean.sa ~ dunif(0, 1) 
mean.p ~ dunif(0, 1)
mean.log.f <- log(mean.f)
mean.f ~ dunif(0, 10)

for (t in 1:(n.occasions-1)){
  p[t] <- mean.p
}

for (t in 1:(n.occasions-1+K)){  # Here we extend the loop to K more years
  logit.sj[t] <- mean.logit.sj + eps.sj[t]
  eps.sj[t] ~ dnorm(0, tau.sj)
  sj[t] <- ilogit(logit.sj[t])
  logit.sa[t] <- mean.logit.sa + eps.sa[t]
  eps.sa[t] ~ dnorm(0, tau.sa)
  sa[t] <- ilogit(logit.sa[t])
}

for (t in 1:(n.occasions+K)){   # Extended loop also here
  log.f[t] <- mean.log.f + eps.f[t]
  eps.f[t] ~ dnorm(0, tau.f)
  f[t] <- exp(log.f[t])
}

sigma.sj ~ dunif(0, 10)
tau.sj <- pow(sigma.sj, -2)
sigma.sa ~ dunif(0, 10)
tau.sa <- pow(sigma.sa, -2)
sigma.f ~ dunif(0, 10)
tau.f <- pow(sigma.f, -2)

sigma ~ dunif(0.5, 50)
tau <- pow(sigma, -2)

# Population count data (state-space model)
# Model for the initial population size: uniform priors
N[1,1] ~ dcat(pNinit)
N[2,1] ~ dcat(pNinit)

# Process model over time: our model of population dynamics
for (t in 1:(n.occasions-1)){  # Note extended loop
  N[1,t+1] ~ dpois(sj[t] * f[t] * (N[1,t] + N[2,t])) 
  N[2,t+1] ~ dbin(sa[t], (N[1,t] + N[2,t]))
}

for (t in n.occasions:(n.occasions-1+K)){  # Note extended loop
  N[1,t+1] ~ dpois(sj.f * f[t] * (N[1,t] + N[2,t])) 
  N[2,t+1] ~ dbin(sa[t], (N[1,t] + N[2,t]))
}

# Observation model
for (t in 1:n.occasions){
  count[t] ~ dnorm(N[1,t] + N[2,t], tau)
}

# Productivity data (Poisson regression model)
for (t in 1:n.occasions){
  J[t] ~ dpois(f[t] * B[t])
}

# Capture-recapture data (CJS model with multinomial likelihood)
# Define the multinomial likelihood
for (t in 1:(n.occasions-1)){
  marr.j[t,1:n.occasions] ~ dmulti(pr.j[t,], rel.j[t])
  marr.a[t,1:n.occasions] ~ dmulti(pr.a[t,], rel.a[t])
}
# Define the cell probabilities of the m-arrays
# Main diagonal
for (t in 1:(n.occasions-1)){
  q[t] <- 1-p[t]   # Probability of non-recapture
  pr.j[t,t] <- sj[t] * p[t]
  pr.a[t,t] <- sa[t] * p[t]
  # Above main diagonal
  for (j in (t+1):(n.occasions-1)){
    pr.j[t,j] <- sj[t] * prod(sa[(t+1):j]) * prod(q[t:(j-1)]) * p[j]
    pr.a[t,j] <- prod(sa[t:j]) * prod(q[t:(j-1)]) * p[j]
  } #j
  # Below main diagonal
  for (j in 1:(t-1)){
    pr.j[t,j] <- 0
    pr.a[t,j] <- 0
  } #j
} #t
# Last column: probability of non-recapture
for (t in 1:(n.occasions-1)){
  pr.j[t,n.occasions] <- 1-sum(pr.j[t,1:(n.occasions-1)])
  pr.a[t,n.occasions] <- 1-sum(pr.a[t,1:(n.occasions-1)])
}

# Derived parameters
# Total population size
for (t in 1:(n.occasions+K)){
  Ntot[t] <- N[1,t] + N[2,t]
}

# Random draw of survival from 1 to n.occassions-1
sj.f <- sample(out1$sims.list$sj, 1)

# Check whether the population is extinct in the future
for (t in 1:K){
  extinct[t] <- equals(Ntot[n.occasions+t], 0)  
}
} # end model
")

# Initial values
inits <- function(){list(mean.sj=runif(1, 0, 0.5), mean.sa=runif(1, 0.4, 0.6), mean.f=runif(1, 1.3, 2))}

# Parameters monitored
parameters <- c("mean.sj", "sigma.sj", "mean.sa", "sigma.sa", "mean.p", "mean.f", "sigma.f", "sigma", "sj", "sa", "f", "N", "Ntot", "extinct")

# MCMC settings
ni <- 20000; nb <- 5000; nc <- 3; nt <- 3; na <- 1000

# Call JAGS (ART 1 min), check convergence and summarize posteriors
out1 <- jags(jags.data, inits, parameters, "model1.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=TRUE)

Answer 1

I think you can randomly draw from the vector of time-dependent survival parameters, but not by sampling from the sims.list object. As @Limey pointed out in the comment above, out1$sims.list is an R object that's not accessible to JAGS and sample is an R function not mentioned in the user manual. Here are the steps:

1. Define a vector of the probabilities for each outcome (probability density function): length n.occasions-1, each with probability 1/(n.occasions-1).
2. Use the categorical distribution, dcat, to randomly sample the density above.
3. Use the random samples as indices to draw samples from the time-dependent survival parameters.

Code:

# Random draw of survival from 1 to n.occassions-1
punif <- rep(1/(n.occasions-1),(n.occasions-1))
b ~ dcat(punif[])
sj.f <- sj[b]

The solution uses this answer, the rep function on page 41 of the JAGS v4.3.0 user manual, and the dcat distribution on page 50 of the user manual (also see this answer).

I could also imagine an alternative solution to this problem. Instead of sampling directly from the estimated time-dependent survival parameters (which are coming from observed years, I think?), you could simulate sj.f from the distribution for survival and introduce a parameter to examine the influence of variability around mean survival. Instead of the lines above, you could test the below, and explore different values for eps.

# Random draw of survival from 1 to n.occassions-1    
logit.sj.f <- mean.logit.sj + eps.sj.f
eps.sj.f ~ dnorm(0, tau.sj+eps)
sj.f <- ilogit(logit.sj.f)