Fred Touzalin
2013-11-06
Dear All,
I am using jags through R2jags on R. I tried to explore the ability of a class of model called "Multievent models" to estimate parameters from simulated data. Here is the error message I got when I try to run jags:
Erreur dans jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
Error in node z[1,7]
Unobserved node inconsistent with unobserved parents at initialization.
Try setting appropriate initial values.
This is very explicit, that is why jags is very friendly, but all my attempts to solve this trouble were unsuccessful...
Before developping my model I began with severals simpler models that works and fit the simulated data I generated. But this one that I thougth straightforward and unambigous seems very sensitive to initail values.
Here is the model:
model { # nind = number of individuals # nocc = number of occasions (observation sessions) #prior phiNBa ~ dunif(0,1) # survival probability for Non-Breeders in site a phiBa ~ dunif(0,1) # survival probability for Breeders in site a phiNBb ~ dunif(0,1) # survival probability for Non-Breeders in site b phiBb ~ dunif(0,1) # survival probability for Breeders in site b psiNBa ~ dunif(0,1) # moving probability for Non-Breeders from site a to site b psiBa ~ dunif(0,1) # moving probability for Breeders from site a to site b psiNBb ~ dunif(0,1) # moving probability for Non-Breeders from site b to site a psiBb ~ dunif(0,1) # moving probability for Breeders from site b to site a Ba ~ dunif(0,1) # breeding probability in site a Bb ~ dunif(0,1) # breeding probability in site b pNBa ~ dunif(0,1) # dectection probability of Non-Breeders in site a pBa ~ dunif(0,1) # dectection probability of Breeders in site a pNBb ~ dunif(0,1) # dectection probability of Non-Breeders in site b pBb ~ dunif(0,1) # dectection probability of Breeders in site b deltaBa ~ dunif(0,1) # uncertainty in the assignment of breeding state in site a deltaBb ~ dunif(0,1) # uncertainty in the assignment of breeding state in site b for (i in 1:nind){ for (t in f[i]:(nocc-1)){ # sate transitions st[1,i,t,1] <- phiNBa*psiNBa*(1-Ba) st[1,i,t,2] <- phiNBa*psiNBa*Ba st[1,i,t,3] <- phiNBa*(1-psiNBa)*(1-Bb) st[1,i,t,4] <- phiNBa*(1-psiNBa)*Bb st[1,i,t,5] <- 1-phiNBa st[2,i,t,1] <- phiBa*psiBa*(1-Ba) st[2,i,t,2] <- phiBa*psiBa*Ba st[2,i,t,3] <- phiBa*(1-psiBa)*(1-Bb) st[2,i,t,4] <- phiBa*(1-psiBa)*Bb st[2,i,t,5] <- 1-phiBa st[3,i,t,1] <- phiNBb*(1-psiNBb)*(1-Ba) st[3,i,t,2] <- phiNBb*(1-psiNBb)*Ba st[3,i,t,3] <- phiNBb*psiNBb*(1-Bb) st[3,i,t,4] <- phiNBb*psiNBb*Bb st[3,i,t,5] <- 1-phiNBb st[4,i,t,1] <- phiBb*(1-psiBb)*(1-Ba) st[4,i,t,2] <- phiBb*(1-psiBb)*Ba st[4,i,t,3] <- phiBb*psiBb*(1-Bb) st[4,i,t,4] <- phiBb*psiBb*Bb st[4,i,t,5] <- 1-phiBb st[5,i,t,1] <- 0 st[5,i,t,2] <- 0 st[5,i,t,3] <- 0 st[5,i,t,4] <- 0 st[5,i,t,5] <- 1 }#t }#i # observation process ob[1,1] <- pNBa ob[1,2] <- 0 ob[1,3] <- 0 ob[1,4] <- 0 ob[1,5] <- 1 - pNBa ob[2,1] <- pBa*(1-deltaBa) ob[2,2] <- pBa*deltaBa ob[2,3] <- 0 ob[2,4] <- 0 ob[2,5] <- 1 - pBa ob[3,1] <- 0 ob[3,2] <- 0 ob[3,3] <- pNBb ob[3,4] <- 0 ob[3,5] <- 1 - pNBb ob[4,1] <- 0 ob[4,2] <- 0 ob[4,3] <- pBb*(1-deltaBb) ob[4,4] <- pBb*deltaBb ob[4,5] <- 1 - pBb ob[5,1] <- 0 ob[5,2] <- 0 ob[5,3] <- 0 ob[5,4] <- 0 ob[5,5] <- 1 # Likelihood for (i in 1:nind){ z[i,f[i]] <- y[i,f[i]] for (t in (f[i]+1):nocc){ z[i,t] ~ dcat(st[z[i,t-1],i,t-1,1:5]) y[i,t] ~ dcat(ob[z[i,t],1:5]) }#t }#i }
I tried to modify the values in the "z" matrix set in jags initial values and I succeeded to make some progress (for expample after several corrections, the error message concerned z[11,5]) until the error message stay the same at the end, whatever the initial values I try. Sometimes I succeeded to set appropriate initial values fo"z", I thought but it was probably not the case, and the error message was no more on "z" but on "y". I also tried to remove "z" matrix in jags data and in jags initial values, but I got a simillar error message.
Maybe the problem is very simple to resolve but I do not succeed to understand it! I am just a jags user and not since so long... and I do not know the way jags check if the initial values are appropriate, that could help.
Many thanks in advance for your suggestions!
Best regards,
Fred
PS: here are the code for the simulating observation data I used ("y" matrix in the code) and to run jags.
####### Simulation ############## # parameters phiNBa = 0.8 phiBa = 0.7 phiNBb = 0.9 phiBb = 0.7 psiNBa = 0.8 psiNBb = 0.7 psiBa = 0.75 psiBb = 0.65 Ba = 0.85 Bb = 0.95 pNBa = 0.8 pBa = 0.6 pNBb= 0.9 pBb = 0.7 deltaBa = 0.2 deltaBb = 0.3 # data creation nind <- 10 nocc <- 12 mst <- matrix(c(phiNBa*psiNBa*(1-Ba), phiNBa*psiNBa*Ba, phiNBa*(1-psiNBa)*(1-Bb), phiNBa*(1-psiNBa)*Bb, 1-phiNBa, phiBa*psiBa*(1-Ba), phiBa*psiBa*Ba, phiBa*(1-psiBa)*(1-Bb), phiBa*(1-psiBa)*Bb, 1-phiBa, phiNBb*(1-psiNBb)*(1-Ba), phiNBb*(1-psiNBb)*Ba, phiNBb*psiNBb*(1-Bb), phiNBb*psiNBb*Bb, 1-phiNBb, phiBb*(1-psiBb)*(1-Ba), phiBb*(1-psiBb)*Ba, phiBb*psiBb*(1-Bb), phiBb*psiBb*Bb, 1-phiBb, 0, 0, 0, 0, 1),nrow=5,byrow=T) round(mst,2); rowSums(mst) mob <- matrix(c(pNBa, 0, 0, 0, 1 - pNBa, pBa*(1-deltaBa), pBa*deltaBa, 0, 0, 1 - pBa, 0, 0, pNBb, 0, 1 - pNBb, 0, 0, pBb*(1-deltaBb), pBb*deltaBb, 1 - pBb, 0, 0, 0, 0, 1), nrow=5,byrow=T) mob; rowSums(mob) data <- matrix(NA,nrow=nind,ncol=nocc) first <- sample(1:5,nind,replace=T) data <- t(data) data[(0:(nind-1))*nocc+first[1:nind]] <- sample(1:4,nind,replace=T) data <- t(data) for (i in 1:nind){ for (t in (first[i]+1):nocc){ data[i,t] <- which(rmultinom(1,1,mst[data[i,t-1],1:5])==1) } } dobs <- matrix(NA,nind,nocc) for (i in 1:nind){ dobs[i,first[i]] <- data[i,first[i]] for (t in (first[i]+1):nocc){ dobs[i,t] <- which(rmultinom(1,1,mob[data[i,t],1:5])==1) } } ####### Jags specifications ############## # data f <- numeric() for (i in 1:dim(dobs)[1]){ f[i] <- min(which(dobs[i,]>=1)) } # known.state known.state <- dobs known.state[known.state %in% c(1,3,5)] <- NA for (i in 1:dim(dobs)[1]){ known.state[i,(1:f[i])] <- NA} known.state #Jags data jags.data <- list(y=dobs, nind=nind, nocc=nocc, f=f, z=known.state) # Initial values # unknown.state unknown.state <- dobs unknown.state[dobs %in% c(2,4)] <- NA x1 <- which(dobs == 1) unknown.state[x1] <- sample(1:2,length(x1),replace=T) x3 <- which(dobs == 3) unknown.state[x3] <- sample(3:4,length(x3),replace=T) x5 <- which(dobs == 5) unknown.state[x5] <- sample(1:5,length(x5),replace=T) for (i in 1:nind){ unknown.state[i,1:f[i]] <- rep(NA,f[i]) } inits <- function(){list(phiNBa=runif(1,0,1),phiBa=runif(1,0,1),phiNBb=runif(1,0,1),phiBb=runif(1,0,1), psiNBa=runif(1,0,1),psiBa=runif(1,0,1),psiNBb=runif(1,0,1),psiBb=runif(1,0,1), Ba=runif(1,0,1),Bb=runif(1,0,1), pNBa=runif(1,0,1),pBa=runif(1,0,1),pNBb=runif(1,0,1),pBb=runif(1,0,1), deltaBa=runif(1,0,1),deltaBb=runif(1,0,1), z = unknown.state)} # Parameters monitored parameters <- c(pNBa, pNBb, pBa, pBb, deltaBa, deltaBb, phiNBa, phiNBb, phiBa, phiBb, psiNBa, psiNBb, psiBa, psiBb,Ba,Bb) # mcmc settings ni <- 1000 nt <- 5 nb <- 500 nc <- 2 # Call JAGS from R library(R2jags) setwd("~/Documents/Avocette/Repro/AnalyseRepro/Simul") start<-as.POSIXlt(Sys.time()) mcmc <- jags(jags.data, inits, parameters, "avo2sites.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, working.directory = getwd()) end <-as.POSIXlt(Sys.time()) duration = end-start; duration
Fred Touzalin
2013-11-07
Hello All,
I am sorry I forgot to say that this model runs with Openbugs without any error message, but I know that Jags is more sensitive to initial values. My preferance is to work with Jags and I would like to understand why jags does not accept the initiales values that I generated !
Best regards,
Fred
PS: there is a mistake in the code at the end of the previous post.
It was written in the section "Parameters monitored" :
parameters <- c(pNBa, pNBb, pBa, pBb, deltaBa, deltaBb, phiNBa, phiNBb, phiBa, phiBb, psiNBa, psiNBb, psiBa, psiBb, Ba, Bb)
but the good code is of course:
parameters <- c("pNBa", "pNBb", "pBa", "pBb", "deltaBa", "deltaBb", "phiNBa", "phiNBb", "phiBa", "phiBb", "psiNBa", "psiNBb", "psiBa", "psiBb", "Ba", "Bb")
Martyn Plummer
2013-11-12
When I run your model I get
Error in node z[1,9] Unobserved node inconsistent with unobserved parents at initialization
Looking at unknown.state
, which contains the initial values for z
, I see
> unknown.state[1,,drop=FALSE] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] NA NA NA 2 2 4 1 5 4 2 4 4
So in the transition from z[1,8] = 5
to z[1,9] = 4
your process escapes from state 5, but your transition matrix st
tells us that 5 is an absorbing state, so this is impossible.
Better to generate initial values that never enter the absorbing state.
Fred Touzalin
2013-11-14
Thank you very much Martyn, it seems that absorbing states should never be use as initial value as you suggest. However in a more complexible model (7 states) the error message after removing the absorbing states from initial values is on y (observation matrix). But this seems due to absorbing states too!
Thank you again for your quick reply and for your very good and usefull program JAGS!