Bad initial values

We have explained how the Gibbs sampler is an algorithm that requires initial values for all the parameters that are being estimated. What often times people do not realize is that, if these initial values are not provided to JAGS, it will generate them by itself, which sometimes can lead to problems.

For example, say we have a simple occupancy model in which \(D_{it}\) is a binary variable indicating if animals of a particular species were seen in location i at time t. Furthermore, let \(z_i\) be a binary variable indicating if location i is occupied or not. We assume that:

\[z_i \sim Bernoulli(\pi)\] \[D_{it} | z_i =1 \sim Bernoulli(\delta)\]

where \(\delta\) is the detection probability and \(\pi\) is the probability of a location being occupied.

In this example, JAGS can run into problems if individuals are observed for a particular location i (i.e., \(D_{it}=1\) ) but JAGS initializes the occupancy status of that location to 0 (i.e., \(z_i=0\)). Based on our modeling assumptions, it would be impossible for us to see an animal if the site is unoccupied and, as a result, JAGS would throw an error! The thing that is most puzzling to people is that sometimes JAGS would run fine and sometimes it would through this error message. The reason for this is that the \(z_i\)’s are being initialized by random draws from the prior, sometimes resulting in sensible \(z_i\)’s and other times resulting in bad \(z_i\)’s. To illustrate this, we simulate some data:

rm(list=ls(all=T))
set.seed(1)

#create some fake data
nloc=100
nsurv=5
pi=0.8
delta=0.5

dat=matrix(0,nloc,nsurv)
z=rbinom(nloc,size=1,prob=pi)
for (i in 1:nloc){
  if (z[i]==1) dat[i,]=rbinom(nsurv,size=1,prob=delta)
}
rownames(dat)=paste('loc',1:nloc)
colnames(dat)=paste('survey',1:nsurv)

#output fake data
# setwd('U:\\uf\\courses\\bayesian course\\rasc\\2016\\chapter15 jags problems\\occupancy example')
setwd('/Volumes/Users/drvalle/uf/courses/bayesian course/rasc/2016/chapter15 jags problems/occupancy example')
write.csv(dat,'fake data occupancy.csv')

Here is the JAGS occupancy model in the file “occupancy model.R”:

model{
  for (i in 1:N){
    Z1[i]~dbern(pi)
    for (j in 1:NSURV){
      D1[i,j]~dbin(ifelse(Z1[i]==0,0,delta),1)  
    }
  }
  pi ~ dbeta(1,1)
  delta ~ dbeta(1,1)
}

Now, you can see this problem happening in the code below:

rm(list=ls(all=T))
library(jagsUI)
set.seed(1)

# MCMC settings 
ni <- 1000 #number of iterations
nt <- 50    #interval to thin 
nb <- ni/2  #number of iterations to discard as burn-in
nc <- 3     #number of chains

#set parameters to monitor
params=c('pi','delta')

#run model
# setwd('U:\\uf\\courses\\bayesian course\\rasc\\2016\\chapter15 jags problems\\occupancy example')
setwd('/Volumes/Users/drvalle/uf/courses/bayesian course/rasc/2016/chapter15 jags problems/occupancy example')
dat=read.csv('fake data occupancy.csv',as.is=T)

model2=jags(model.file="occupancy model.R",
            parameters.to.save=params,data=list(D1=dat[,-1],N=nrow(dat),NSURV=ncol(dat)-1),
            n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE)
## 
## Processing function input....... 
## 
## Converted data.frame D1 to matrix
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 500
##    Unobserved stochastic nodes: 102
##    Total graph size: 806
## 
## Initializing model
## Deleting model
## Error in rjags::jags.model(file = model.file, data = data, inits = inits, : Error in node D1[1,1]
## Node inconsistent with parents

To solve this problem, we have to initialize the \(z_i\)’s with better values. For instance, we can set \(z_i=1\) for all locations. This is fine even if the corresponding \(D_{it}=0\) for all t in that location. Alternatively, we can set \(z_i=1\) for all locations in which at least one \(D_{it}\) was equal to 1. Notice that we are just setting the initial values for these random variables and therefore these initial values should (in principle) be inconsequential. The code below illustrates how this solution solves the problem we experienced before. I just modify a couple of lines in my code to make this change (emphasized with comment “NEW”):

rm(list=ls(all=T))
library(jagsUI)
set.seed(1)

# MCMC settings 
ni <- 1000 #number of iterations
nt <- 50    #interval to thin 
nb <- ni/2  #number of iterations to discard as burn-in
nc <- 3     #number of chains

#set parameters to monitor
params=c('pi','delta')

#run model
# setwd('U:\\uf\\courses\\bayesian course\\rasc\\2016\\chapter15 jags problems\\occupancy example')
setwd('/Volumes/Users/drvalle/uf/courses/bayesian course/rasc/2016/chapter15 jags problems/occupancy example')
dat=read.csv('fake data occupancy.csv',as.is=T)

#NEW
inits<-function(){
  list(Z1=rep(1,nrow(dat)))
}

model2=jags(model.file="occupancy model.R",
            inits=inits,                    #NEW
            parameters.to.save=params,data=list(D1=dat[,-1],N=nrow(dat),NSURV=ncol(dat)-1),
            n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE)
## 
## Processing function input....... 
## 
## Converted data.frame D1 to matrix
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 500
##    Unobserved stochastic nodes: 102
##    Total graph size: 806
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 500 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 500 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
#check to see if results make sense
summary(model2)
##                 mean          sd        2.5%         25%         50%
## pi         0.8367588  0.04139626   0.7638299   0.8133039   0.8358745
## delta      0.4604549  0.02545876   0.4123724   0.4490819   0.4635433
## deviance 584.0553997 11.33012511 561.0902080 579.6419020 585.8126895
##                  75%       97.5%      Rhat n.eff overlap0 f
## pi         0.8652467   0.9091295 0.9607404    30        0 1
## delta      0.4777580   0.4981036 0.9763118    30        0 1
## deviance 592.2564231 600.6382714 1.0198678    30        0 1



Comments?

Send me an email at

References