The R package jagsUI is a wrapper for JAGS that has some awesome functions, including a posterior predictive check. As discussed here, you simulate the new data for the parameters based on the posterior distribution and then compare these to the actual data to determine how well the model fits. https://rdrr.io/cran/jagsUI/man/ppcheck.html#heading-3

However, I am having trouble expanding this to my more complex hierarchical model. I have a multi-state Cormack-Jolly-Seber model using capture recapture data that estimates survival and recapture probability of juveniles and adults in my system at two different locations, high and low.

My data is in a matrix format of x=528 birds y=4 years (also referred to as n.occasions), in which 1= seen as juvenile, 2= seen as adult, and 3= unseen. So a bird that was seen as a juvenile in the first year and never seen again would have a capture history of 1333 while a bird seen in the first year as an adult and seen only the following year would be 2233. A juvenile first seen in the second year of the study and sighted again as an adult the next 2 years would be 3122.

Site is included as a fixed group effect (called "group") by including a vector with length = # birds in which 1= high site and 2= low site.

Please does anybody have advice on how to change my model to incorporate a ppcheck?

Code originally sourced from Bayesian Population Analysis Using WinBUGS by Kery and Schaub 2012.

Vector with occasion of first capture

get.first <- function(x) min(which(x!=0)) f <- apply(data,1,get.first)

Recode data matrix so that 3 = not seen

data[data==0] <- 3

Analysis of the model

sink("ageg.jags") cat(" model { #---------------------------------------- # Parameters: # phi.juv: juvenile survival probability # phi.ad: adult survival probability # p: recapture probability #---------------------------------------- # States (S): # 1 alive as juvenile # 2 alive as adult # 3 dead # Observations (O): # 1 seen as juvenile # 2 seen as adult # 3 not seen #----------------------------------------

# Priors and constraints
for(i in 1:nind){
for (t in f[i]:(n.occasions-1)){
phi.juv[i,t] <- mean.phijuv.g[group[i]]
phi.ad[i,t] <- mean.phiad.g[group[i]]
p[i,t] <- mean.p.g[group[i]]
} #t
for (u in 1:g){
mean.phijuv.g[u] ~ dunif(0, 1)     #priors for group-specific juvenile survival
mean.phiad.g[u] ~ dunif(0, 1)      #priors for group-specific adult survival
mean.p.g[u] ~ dunif(0, 1)          #priors for group-specific recapture

# Define state-transition and observation matrices  
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in f[i]:(n.occasions-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phi.juv[i,t]
ps[1,i,t,3] <- 1-phi.juv[i,t]
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phi.ad[i,t]
ps[2,i,t,3] <- 1-phi.ad[i,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] <- 0
po[1,i,t,3] <- 1
po[2,i,t,1] <- 0
po[2,i,t,2] <- p[i,t]
po[2,i,t,3] <- 1-p[i,t]
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- 1
} #t
} #i

# State-space model likelihood
for (i in 1:nind){
z[i,f[i]] <- y[i,f[i]]
for (t in (f[i]+1):n.occasions){
# State equation: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1,])
# Observation equation: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1,])
} #t
} #i
} #model
",fill = TRUE)


Bundle data

jags.data <- list(y = data, f = f, n.occasions = dim(data)[2], nind = dim(data)[1],g=length(unique(group)), group=group)

Initial values function

age.ms.init <- function(ch){ init <- ch al <- which(ch==2, arr.ind = T) for (i in 1:nrow(al)){ init[al[i,1], f[al[i,1]]:al[i,2]] <- 2 } #i for (i in 1:nrow(init)){ init[i,1:f[i]] <- NA } #i return(init) }

initial values

inits <- function(){list(mean.phijuv.g = runif(length(unique(group)), 0, 1), mean.phiad.g = runif(length(unique(group)), 0, 1), mean.p.g = runif(length(unique(group)), 0, 1), z = age.ms.init(data))}

Parameters monitored, phi is for survival, p is for recapture probability

parameters <- c("mean.phijuv.g", "mean.phiad.g", "mean.p.g")

MCMC settings

ni <- 2000 nt <- 3 nb <- 1000 nc <- 3 nad <- 1000

Call JAGS using jagsUI

age_g <- jags(jags.data, inits, parameters, "ageg.jags",n.chains = nc, n.adapt = nad, n.thin = nt, n.iter = ni, n.burnin = nb)

The output is 4 different apparent survival estimates (high juvenile, high adult, low juvenile, low adult) and then estimates of recapture probability for both sites (high and low).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.