我正在尝试使用 Line-Transect Distance Sampling 和数据增强来模拟乌龟洞穴的大小。但是,我不断收到错误消息“无法找到合适的采样器”。
一些背景:洞穴的宽度可以从 4 厘米到 55 厘米不等,并且以各种概率和与景观中的观察者不同的距离被看到。为了实现收敛,我决定使用基于分类的模型。乌龟洞穴被放入 7 个箱中的一个,在给定箱中的概率来自 Dirichlet 分布。
在这部分代码中,增强的洞穴从模型中提取它们的大小:
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
}
一旦它们有了大小,就确定了看到洞穴 p[i] 的概率:
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
所有这些似乎都可以正常工作,但是在评估 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])
我可以运行除 y 之外的每一行,并且模型运行良好。但是添加该行并包含我的 y 数据会引发“无法找到合适的采样器”错误。知道为什么吗?我需要包含这个变量,因为 y 表示我是否真的在调查期间找到了洞穴。
任何建议将不胜感激。
这是带有一些虚假数据的完整代码:
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
}
"
数据:
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))
运行命令:
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)