Introduction

One thing is to listen to these problems and hear how one can solve them. Another very different thing is actually solving problems in the real life.

Problem 1

Imagine that we are interested in estimating the true number of individuals at a particular location \(N_i\) and that we have data on the number of individuals that we saw at location i in survey t \(X_{it}\). Further assume that the same individual cannot be counted multiple times within survey t and that population size does not change through time (this assumption should be fine if the surveys are done in a short period of time). In this case, perhaps a reasonable model would be:

\[X_{it} \sim Binomial(N_i,\delta)\]

where \(\delta\) is the detection probability, assumed to be the same for all individuals, all locations, and all surveys. Furthermore, say that we want to understand how population size \(N_i\) is related to vegetation. To this end, we add another layer to our model by assuming a Poisson regression:

\[N_i \sim Poisson(exp(\beta_0 + \beta_1 veg_i))\]

This is a N-mixture model, a type of model widely used in the context of wildlife studies. We define our model in the file “pop size model.R”:

model{
  for (i in 1:nloc){
    N1[i] ~ dpois(exp(b0+b1*veg[i]))
    for (j in 1:nsurv){
      D1[i,j]~dbin(N1[i],delta)  
    }
  }
  delta ~ dbeta(1,1)
  b0 ~ dnorm(0,0.1)
  b1 ~ dnorm(0,0.1)
}

The data to fit this model is called “fake data pop size.csv”. We attempt to fit this model using the following code:

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('N1','delta','b0','b1')

#run model
# setwd('/Volumes/Users/drvalle/uf/courses/bayesian course/group activities/20 debug')
setwd('U:\\uf\\courses\\bayesian course\\group activities\\20 debug')
dat=read.csv('fake data pop size.csv',as.is=T)
ind=grep('survey',colnames(dat))
D1=data.matrix(dat[,ind])
veg=dat$veg

#set initial values
inits=function(){
  list(N1=rpois(nrow(D1),lambda=5),  
       delta=runif(1),
       b0=rnorm(1,0,10),
       b1=rnorm(1,0,10))
}

#fit model
model2=jags(model.file="pop size model.R",
            parameters.to.save=params,
            data=list(D1=D1,nloc=nrow(D1),
                      nsurv=ncol(D1),veg=veg),
            n.chains=nc,n.burnin=nb,n.iter=ni,n.thin=nt,DIC=TRUE)

model2

How come JAGS runs into problems with this model? Once you have identified the problem(s), please fix it to make sure that everything is working properly.

Back to main menu