Initial values for hidden process in a multinomial model

Help
2013-11-06
2013-11-14
  • Fred Touzalin
    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
    
     
    Last edit: Martyn Plummer 2013-11-12
  • Fred Touzalin
    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
    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
    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!