I'm trying to model the sizes of tortoise burrows using Line-Transect Distance Sampling and data augmentation. However, I keep getting the error "unable to find appropriate sampler."
Some background: Burrows can be from 4cm to 55cm wide and are seen with various probabilities and at various distances from the observer in the landscape. In order to achieve convergence, I've decided to go with a categorical based model. Tortoise burrows are put into one of 7 bins, with the probability of being in a given bin coming from a Dirichlet distribution.
Augmented burrows draw their sizes from the model in this part of the code:
for (i in 1:(nind +nz)) {
z[i] ~ dunif(a[bin.size[i]],b[bin.size[i]]) #augmented burrows need real sizes; drawn from bin constraints
bins[i, ] ~ dmulti(pClust[1:7],1)
bin.size[i] <- sum(bins[i,1]+bins[i,2]*2+bins[i,3]*3+bins[i,4]*4+bins[i,5]*5+bins[i,6]*6+bins[i,7]*7)
}
pClust[1:7] ~ ddirch(psizes) #probability of being in a bin
for (ii in 1:7) {
psizes[ii]~ dunif(0,1)
}
for (clustIdx in 1:1) {
a[clustIdx] <- 4 #the smallest possible burrows are 4cm; smaller than that makes no biological sense
b[clustIdx] <- 9.9
}
for (clustIdx in 2:2) {
a[clustIdx] <- 10
b[clustIdx] <- 17.9
}
for (clustIdx in 3:3) {
a[clustIdx] <- 18
b[clustIdx] <- 22.9
}
for (clustIdx in 4:4) {
a[clustIdx] <- 23
b[clustIdx] <- 28.9
}
for (clustIdx in 5:5) {
a[clustIdx] <- 29
b[clustIdx] <- 34.9
}
for (clustIdx in 6:6) {
a[clustIdx] <- 35
b[clustIdx] <- 40.9
}
for (clustIdx in 7:7) {
a[clustIdx] <- 41
b[clustIdx] <- 55
}
Once they have a size, the probability of seeing that burrow p[i] is determined:
sigma[i] <- exp(sigma.int+sigma.beta*z[i]) #log link for size
logp[i] <- -((x[i]*x[i])/(2*sigma[i]*sigma[i])) #This is the normal distribution with an adjustment for covs
p[i] <- exp(logp[i])*xi[i] # probabilty of seeing the thing, regardless of us actually seeing it; scaled by xi
xi[i] <- ifelse(z[i] < b.point, m*z[i]+intercept, 1) #if less than b.point, probability is linear; if larger than b.point, perfect detection on line
}
m <- (1-p.online)/(b.point-4) #slope for detection on the line for smaller burrows
intercept <- p.online-(4*m) ## finding intercept via the detection of the 4cm burrow
p.online ~dunif(0.2, 0.5) #estimated detection on line for 4 cm burrows
b.point~dunif(10,30) #where the stick breaks
All of that seems to work fine, but the problem happens here when evaluating y:
mu[i] <- w[i]*p[i] ## probabilty of seeing it for all the ones we DID see (where w[i] = 1)
y[i] ~ dbern(mu[i])
I can run every line excluding the y's and the model runs fine. But adding in that line and including my y data throws the "unable to find appropriate sampler" error. Any idea why? I need this variable to be included because y represents if I actually found the burrow during our survey or not.
Any suggestions would be greatly appreciated.
Here's the full code with some fake data:
modelstring.NoIntense = "
model
{
for (i in 1:(nind +nz)) {
z[i] ~ dunif(a[bin.size[i]],b[bin.size[i]]) #augmented burrows need real sizes; drawn from bin constraints
bins[i, ] ~ dmulti(pClust[1:7],1)
bin.size[i] <- sum(bins[i,1]+bins[i,2]*2+bins[i,3]*3+bins[i,4]*4+bins[i,5]*5+bins[i,6]*6+bins[i,7]*7)
w[i] ~ dbern(psi) ##augmentation
x[i] ~ dunif(0,Bx) #distance from line for the missed ones; Bx = max(distances)
sigma[i] <- exp(sigma.int+sigma.beta*z[i]) #log link for size
logp[i] <- -((x[i]*x[i])/(2*sigma[i]*sigma[i])) #This is the normal distribution with an adjustment for covs
p[i] <- exp(logp[i])*xi[i] # probabilty of seeing the thing, regardless of us actually seeing it; scaled by xi
xi[i] <- ifelse(z[i] < b.point, m*z[i]+intercept, 1) #if less than b.point, probability is linear; if larger than b.point, perfect detection on line
mu[i] <- w[i]*p[i] ## probabilty of seeing it for all the ones we DID see (where w[i] = 1)
y[i] ~ dbern(mu[i])
}
m <- (1-p.online)/(b.point-4) #slope for detection on the line for smaller burrows
intercept <- p.online-(4*m) ## finding intercept via the detection of the 4cm burrow
p.online ~dunif(0.2, 0.5) #estimated detection on line for 4 cm burrows
b.point~dunif(10,30) #where the stick breaks
pClust[1:7] ~ ddirch(psizes) #probability of being in a bin
for (ii in 1:7) {
psizes[ii]~ dunif(0,1)
}
for (clustIdx in 1:1) {
a[clustIdx] <- 4 #the smallest possible burrows are 4cm; smaller than that makes no biological sense
b[clustIdx] <- 9.9
}
for (clustIdx in 2:2) {
a[clustIdx] <- 10
b[clustIdx] <- 17.9
}
for (clustIdx in 3:3) {
a[clustIdx] <- 18
b[clustIdx] <- 22.9
}
for (clustIdx in 4:4) {
a[clustIdx] <- 23
b[clustIdx] <- 28.9
}
for (clustIdx in 5:5) {
a[clustIdx] <- 29
b[clustIdx] <- 34.9
}
for (clustIdx in 6:6) {
a[clustIdx] <- 35
b[clustIdx] <- 40.9
}
for (clustIdx in 7:7) {
a[clustIdx] <- 41
b[clustIdx] <- 55
}
sigma.int~ dnorm(0,s.int.tau)
s.int.tau <- 1/(s.int.sd*s.int.sd)
s.int.sd ~ dunif(.00001,5)
sigma.beta~ dnorm(0,s.beta.tau)T(0,)
s.beta.tau <- 1/(s.beta.sd*s.beta.sd)
s.beta.sd ~ dunif(.00001,5)
o.int~ dnorm(0,o.int.tau)
o.int.tau <- 1/(o.int.sd*o.int.sd)
o.int.sd ~ dunif(.00001,5)
psi~ dunif(0,1) #exists or not
N <- sum(w)
D <- N/(2*L*Bx) #burrow density
}
"
Data:
nind <- 100
nz <- 400
z <- c(rnorm(70,32, 5), rnorm(10, 10, 3), rnorm(20,40,2), rep(NA,400))
z2 <- matrix(NA, ncol = 7, nrow = nind+nz)
for (kk in 1:(nind+nz)){
z2[kk,] <- ifelse(rep(z[kk] < 10,7), (c(1,rep(0,6))),
ifelse(rep(z[kk] <18 & z[kk] >= 10,7), (c(0,1,rep(0,5))),
ifelse(rep(z[kk] >= 18 & z[kk] <23,7), (c(0,0,1, rep(0,4))),
ifelse(rep(z[kk] >=23 & z[kk] <29,7), (c(rep(0,3),1,0,0,0)),
ifelse(rep(z[kk] >=29& z[kk] <35,7), (c(rep(0,4), 1, 0, 0)),
ifelse(rep(z[kk] >=35 & z[kk] <41,7), (c(rep(0,5), 1, 0)),
ifelse(rep(z[kk] >= 41,7), (c(rep(0,6),1)), NA
)))))))
}
psizes = vector()
for (i in 1:7) {
psizes[i] = sum(z2[,i], na.rm = T)
}
L = 1.2
x = c(abs(rnorm(70,0,10)), abs(rnorm(10,0,4)), abs(rnorm(20,0,15)), rep(NA,400))
Bx = max(x, na.rm = TRUE)
y = c(rep(1,100), rep(0,400))
o = c(rbinom(100, 1, .75), rep (NA,400))
w = c(rep(1,100), rep(NA,400))
Run Commands:
jd.NoIntense = list(nind= nind, nz = nz, z= z, bins = z2, L = L, Bx = Bx, x=x,o=o, w= w,y=y, psizes = psizes)
ji.NoIntense <- list(pClust = psizes/nind,sigma.int = -.02, sigma.beta = 2, o.int = .01, p.online = .35, b.point = 20)
jp.NoIntense <- c("pClust","sigma.int", "sigma.beta", "p.online", "b.point", "o.int", "p[101]", "N")
NoIntense <- run.jags(modelstring.NoIntense, data = jd.NoIntense, inits = list(ji.NoIntense,ji.NoIntense, ji.NoIntense), monitor= jp.NoIntense,
burnin = 10000, sample = 5000, adapt = 1000, method = 'parallel', n.chain = 3, psrf.target = 1.1)