Я пытался запустить многоуровневую модель метки-повторного захвата, модифицированную из примеров в Bayesian Population Analysis, используя WinBUGS от Kery and Schaub (2012). Тем не менее, я застрял на некоторое время, поэтому я надеюсь, что кто-то может иметь некоторое понимание. Я продолжаю получать следующую ошибку:
Error in jags.model(file = model.file, data = data, inits = inits, n.chains = n.chains, :
Error in node y[1,7]
Node inconsistent with parents
Я довольно новичок в JAGS (и в байесовском анализе в целом), но из того, что я прочитал, я понимаю, что это не редкая ошибка и проблема скорее всего лежит с начальными значениями, которые я предоставил. Я пробовал различные начальные значения, но продолжаю получать ту же ошибку. Я также пробовал гамма-априоры для вероятностей перехода состояний (psi), но получаю ту же ошибку. Моя переменная ответа (y) представляет собой матрицу 392 x 99, что означает 99 случаев выборки, что довольно мало - интересно, связана ли проблема с количеством случаев выборки? (Когда я запускаю пример точно так же, как в Kery and Schaub (2012) с 6 случаями выборки, он работает.) Это симулированные данные, но я получаю то же сообщение об ошибке с моими реальными данными, которые расположены в матрице аналогичного размера. Я предоставил код для имитации данных и код моей модели ниже. Обратите внимание, что вероятность перехода из состояний BF в NB установлена равной нулю - индивидуумы не переходят в этом направлении. Если у кого-то есть идеи относительно более подходящих начальных значений или другого источника проблемы, я был бы очень признателен.
################################################################################
# MULTISTATE MARK-RECAPTURE MODEL
# adapted from Kery and Schaub (2012)
################################################################################
#### SIMULATE DATA ####
library(rjags)
library(jagsUI)
set.seed(42)
# Define mean survival, transitions, recapture, as well as number of occasions, states, observations and released individuals
phiNB <- 0.85
phiNST <- 0.75
phiBF <- 0.65
psiNB.NST <- 0.3
psiNB.BF <- 0.2
psiNST.NB <- 0.5
psiNST.BF <- 0.1
psiBF.NB <- 0
psiBF.NST <- 0.3
psiNB <- 0.7
psiNST <- 0.4
psiBF <- 0.5
n.occasions <- 99
n.states <- 4
n.obs <- 4
marked <- matrix(NA, ncol = n.states, nrow = n.occasions)
marked[,1] <- rep(1, n.occasions)
marked[,2] <- rep(2, n.occasions)
marked[,3] <- rep(1, n.occasions)
marked[,4] <- rep(0, n.occasions)
# Define matrices with survival, transition and recapture probabilities
# These are 4-dimensional matrices, with
# Dimension 1: state of departure
# Dimension 2: state of arrival
# Dimension 3: individual
# Dimension 4: time
# 1. State process matrix
totrel <- sum(marked)*(n.occasions-1)
PSI.STATE <- array(NA, dim=c(n.states, n.states, totrel, n.occasions-1))
for (i in 1:totrel){
for (t in 1:(n.occasions-1)){
PSI.STATE[,,i,t] <- matrix(c(
phiNB*(1-psiNB.NST-psiNB.BF), phiNB*psiNB.NST, phiNB*psiNB.BF, 1-phiNB,
phiNST*psiNST.NB, phiNST*(1-psiNST.NB-psiNST.BF), phiNST*psiNST.BF, 1-phiNST,
phiBF*psiBF.NB, phiBF*psiBF.NST, phiBF*(1-psiBF.NB-psiBF.NST), 1-phiBF,
0, 0, 0, 1), nrow = n.states, byrow = TRUE)
} #t
} #i
# 2.Observation process matrix
PSI.OBS <- array(NA, dim=c(n.states, n.obs, totrel, n.occasions-1))
for (i in 1:totrel){
for (t in 1:(n.occasions-1)){
PSI.OBS[,,i,t] <- matrix(c(
psiNB, 0, 0, 1-psiNB,
0, psiNST, 0, 1-psiNST,
0, 0, psiBF, 1-psiBF,
0, 0, 0, 1), nrow = n.states, byrow = TRUE)
} #t
} #i
# Define function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, marked, unobservable = NA){
# Unobservable: number of state that is unobservable
n.occasions <- dim(PSI.STATE)[4] + 1
CH <- CH.TRUE <- matrix(NA, ncol = n.occasions, nrow = sum(marked))
# Define a vector with the occasion of marking
mark.occ <- matrix(0, ncol = dim(PSI.STATE)[1], nrow = sum(marked))
g <- colSums(marked)
for (s in 1:dim(PSI.STATE)[1]){
if (g[s]==0) next # To avoid error message if nothing to replace
mark.occ[(cumsum(g[1:s])-g[s]+1)[s]:cumsum(g[1:s])[s],s] <-
rep(1:n.occasions, marked[1:n.occasions,s])
} #s
for (i in 1:sum(marked)){
for (s in 1:dim(PSI.STATE)[1]){
if (mark.occ[i,s]==0) next
first <- mark.occ[i,s]
CH[i,first] <- s
CH.TRUE[i,first] <- s
} #s
for (t in (first+1):n.occasions){
# Multinomial trials for state transitions
if (first==n.occasions) next
state <- which(rmultinom(1, 1, PSI.STATE[CH.TRUE[i,t-1],,i,t-1])==1)
CH.TRUE[i,t] <- state
# Multinomial trials for observation process
event <- which(rmultinom(1, 1, PSI.OBS[CH.TRUE[i,t],,i,t-1])==1)
CH[i,t] <- event
} #t
} #i
# Replace the NA and the highest state number (dead) in the file by 0
CH[is.na(CH)] <- 0
CH[CH==dim(PSI.STATE)[1]] <- 0
CH[CH==unobservable] <- 0
id <- numeric(0)
for (i in 1:dim(CH)[1]){
z <- min(which(CH[i,]!=0))
ifelse(z==dim(CH)[2], id <- c(id,i), id <- c(id))
}
return(list(CH=CH[-id,], CH.TRUE=CH.TRUE[-id,]))
# CH: capture histories to be used
# CH.TRUE: capture histories with perfect observation
}
# Execute simulation function
sim <- simul.ms(PSI.STATE, PSI.OBS, marked)
CH <- sim$CH
# Compute vector with occasions of first capture
get.first <- function(x) min(which(x!=0))
f <- apply(CH, 1, get.first)
# Recode CH matrix: note, a 0 is not allowed in WinBUGS! (is it in JAGS?)
# 1 = seen alive in NB, 2 = seen alive in NST, 3, seen alive in BF, 4 = not seen
rCH <- CH # Recoded CH
rCH[rCH==0] <- 4
#### FUNCTIONS FOR INITS ####
# Function to create known latent states z
known.state.ms <- function(ms, notseen){
# notseen: label for ‘not seen’
state <- ms
state[state==notseen] <- NA
for (i in 1:dim(ms)[1]){
m <- min(which(!is.na(state[i,])))
state[i,m] <- NA
}
return(state)
}
# Function to create initial values for unknown z
ms.init.z <- function(ch, f){
for (i in 1:dim(ch)[1]){ch[i,1:f[i]] <- NA}
states <- max(ch, na.rm = TRUE)
known.states <- 1:(states-1)
v <- which(ch==states)
ch[-v] <- NA
ch[v] <- sample(known.states, length(v), replace = TRUE)
return(ch)
}
#### JAGS MODEL CODE ####
sink("SESA-multistate.jags")
cat("
model {
# -------------------------------------------------
# Parameters:
# phiNB: persistence (survival) probability at site NB
# phiNST: persistence (survival) probability at site NST
# phiBF: persistence (survival) probability at site BF
# psiNB.NST: movement (transition) probability from NB to NST
# psiNB.BF: movement (transition) probability from NB to BF
# psiNST.NB: movement (transition) probability from NST to NB
# psiNST.BF: movement (transition) probability from NST to BF
# psiBF.NB: movement (transition) probability from BF to NB
# psiBF.NST: movement (transition) probability from BF to NST
# pNB: detection probability at site NB
# pNST: detection probability at site NST
# pBF: detection probability at site BF
# -------------------------------------------------
# States (S):
# 1 alive at site NB
# 2 alive at site NST
# 3 alive at site BF
# 0 gone
# Observations (O):
# 1 detected at site NB
# 2 detected at site NST
# 3 detected at site BF
# 0 not detected
# -------------------------------------------------
# Priors and constraints
# Survival and recapture: uniform
phiNB ~ dunif(0, 1)
phiNST ~ dunif(0, 1)
phiBF ~ dunif(0, 1)
pNB ~ dunif(0, 1)
pNST ~ dunif(0, 1)
pBF ~ dunif(0, 1)
# Transitions: multinomial logit
# Normal priors on logit of all but one transition probas
for (i in 1:2){
lpsiNB[i] ~ dnorm(0, 0.001)
lpsiNST[i] ~ dnorm(0, 0.001)
lpsiBF[i] ~ dnorm(0, 0.001)
}
# Constrain the transitions such that their sum is < 1
for (i in 1:2){
psiNB[i] <- exp(lpsiNB[i]) / (1 + exp(lpsiNB[1]) + exp(lpsiNB[2]))
psiNST[i] <- exp(lpsiNST[i]) / (1 + exp(lpsiNST[1]) + exp(lpsiNST[2]))
psiBF[i] <- exp(lpsiBF[i]) / (1 + exp(lpsiBF[1]) + exp(lpsiBF[2]))
}
# Calculate the last transition probability
psiNB[3] <- 1-psiNB[1]-psiNB[2]
psiNST[3] <- 1-psiNST[1]-psiNST[2]
psiBF[3] <- 1-psiBF[1]-psiBF[2]
# 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] <- phiNB * psiNB[1]
ps[1,i,t,2] <- phiNB * psiNB[2]
ps[1,i,t,3] <- phiNB * psiNB[3]
ps[1,i,t,4] <- 1-phiNB
ps[2,i,t,1] <- phiNST * psiNST[1]
ps[2,i,t,2] <- phiNST * psiNST[2]
ps[2,i,t,3] <- phiNST * psiNST[3]
ps[2,i,t,4] <- 1-phiNST
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiBF * psiBF[2]
ps[3,i,t,3] <- phiBF * psiBF[3]
ps[3,i,t,4] <- 1-phiBF
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- pNB
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 1-pNB
po[2,i,t,1] <- pNST
po[2,i,t,2] <- 0
po[2,i,t,3] <- 0
po[2,i,t,4] <- 1-pNST
po[3,i,t,1] <- pBF
po[3,i,t,2] <- 0
po[3,i,t,3] <- 0
po[3,i,t,4] <- 1-pBF
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,f[i]] <- y[i,f[i]]
for (t in (f[i]+1):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
}
",fill = TRUE)
sink()
#### R CODE TO RUN MODEL ####
# Bundle data
jags.data <- list(y = rCH, f = f, n.occasions = dim(rCH)[2], nind = dim(rCH)[1],
z = known.state.ms(rCH, 4))
# Initial values
inits <- function(){list(phiNB = runif(1, 0, 1), phiNST = runif(1, 0, 1), phiBF = runif(1, 0, 1),
lpsiNB = rnorm(2), lpsiNST = rnorm(2), lpsiBF = rnorm(2),
pNB = runif(1, 0, 1) , pNST = runif(1, 0, 1) , pBF = runif(1, 0, 1),
z = ms.init.z(rCH, f))}
# Parameters monitored
parameters <- c("phiNB", "phiNST", "phiBF",
"psiNB", "psiNST", "psiBF",
"pNB", "pNST", "pBF")
# MCMC settings
ni <- 100 # very few just to see if model will initialize
nt <- 1
nb <- 20
nc <- 3
# Call JAGS from R (BRT 56 min)
ms.test <- jags(jags.data, inits, parameters, "SESA-multistate.jags",
n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
print(ms.test, digits = 3)
Заранее большое спасибо!
(Примечание: я я публикую это на справочном форуме JAGS )