我是 Royle-Nichols 多类群占用模型的新手。运行我在此处附加的代码后,我得到“jags.model 中的错误(model.file,data = data,inits = init.values,n.chains = n.chains,:节点 y[2,1] 节点中的错误与父母不合”。
我在此处发布的相关问题中看到该问题可能与我的先验和/或缩写的大小有关...我尝试了一些更改,但问题仍然存在。
如果有人可以在这里提供帮助,我将不胜感激。
sink(file = here("data","RN.model.formula.aug.txt"))
cat("model {
# prior distributions on community level estimates - hyperparameters
psi ~ dunif(0,1) # inclusion rate that generates wi
# mean value (mu)
# parameter related to abundance
mu.a0 ~ dnorm(0,0.5) # intercept on lambda
mu.a1 ~ dnorm(0,0.5) # slope on lambda for hydrography
mu.a2 ~ dnorm(0,0.5) # slope on lambda for floodded
mu.a3 ~ dnorm(0,0.5) # slope on lambda for deforestation
mu.a4 ~ dnorm(0,0.5) # slope on lambda for PA
mu.a5 ~ dnorm(0,0.5) # slope on lambda for City HP
mu.a6 ~ dnorm(0,0.5) # slope on lambda for Local HP
mu.com~ dnorm(0,0.5) # intecept on lambda for community
# parameter related to detectability
mu.r0 ~ dnorm(0,0.5) # intercept on lambda
mu.r1 ~ dnorm(0,0.5) # slope on lambda for hydrography
mu.r2 ~ dnorm(0,0.5) # slope on lambda for floodded
mu.r3 ~ dnorm(0,0.5) # slope on lambda for deforestation
mu.r4 ~ dnorm(0,0.5) # slope on lambda for PA
mu.r5 ~ dnorm(0,0.5) # slope on lambda for City HP
mu.r6 ~ dnorm(0,0.5) # slope on lambda for Local HP
mu.r7 ~ dnorm(0,0.5) # slope on lambda for Effort
# standard deviation
# parameter related to abundance
sigma.a0 ~ dunif(0,10) # intercept
sigma.a1 ~ dunif(0,10) # hydrography
sigma.a2 ~ dunif(0,10) # floodded
sigma.a3 ~ dunif(0,10) # deforestation
sigma.a4 ~ dunif(0,10) # PA
sigma.a5 ~ dunif(0,10) # City HP
sigma.a6 ~ dunif(0,10) # Local HP
sigma.com~ dunif(0,10) # Community
# parameter related to detectability
sigma.r0 ~ dunif(0,10) # intercept
sigma.r1 ~ dunif(0,10) # hydrography
sigma.r2 ~ dunif(0,10) # floodded
sigma.r3 ~ dunif(0,10) # deforestation
sigma.r4 ~ dunif(0,10) # PA
sigma.r5 ~ dunif(0,10) # City HP
sigma.r6 ~ dunif(0,10) # Local HP
sigma.r7 ~ dunif(0,10) # Effort
# create precision
# parameter related to abundance
tau.a0 <- pow(sigma.a0,-2)
tau.a1 <- pow(sigma.a1,-2)
tau.a2 <- pow(sigma.a2,-2)
tau.a3 <- pow(sigma.a3,-2)
tau.a4 <- pow(sigma.a4,-2)
tau.a5 <- pow(sigma.a5,-2)
tau.a6 <- pow(sigma.a6,-2)
tau.com<- pow(sigma.com,-2)
# parameter related to detectability
tau.r0 <- pow(sigma.r0,-2)
tau.r1 <- pow(sigma.r1,-2)
tau.r2 <- pow(sigma.r2,-2)
tau.r3 <- pow(sigma.r3,-2)
tau.r4 <- pow(sigma.r4,-2)
tau.r5 <- pow(sigma.r5,-2)
tau.r6 <- pow(sigma.r6,-2)
tau.r7 <- pow(sigma.r7,-2)
## generating priors for Community random variable for each species; governed by community-level hyperparameters
for(i in 1:(nspecies+nzeros)){
for (b in 1:n.comunidade){
acom[b,i] ~ dnorm (mu.com, tau.com) # random effect communities
} #j
} #i
for(i in 1:(nspecies+nzeros)) {
# generating parameters for each species related to abundance; governed by community-level hyperparameters
a0[i] ~ dnorm(mu.a0,tau.a0)#I(-10,10)
a1[i] ~ dnorm(mu.a1,tau.a1)#I(-10,10)
a2[i] ~ dnorm(mu.a2,tau.a2)#I(-10,10)
a3[i] ~ dnorm(mu.a3,tau.a3)#I(-10,10)
a4[i] ~ dnorm(mu.a4,tau.a4)#I(-10,10)
a5[i] ~ dnorm(mu.a5,tau.a5)#I(-10,10)
a6[i] ~ dnorm(mu.a6,tau.a6)#I(-10,10)
# generating parameters for each species related to individual detection; governed by community-level hyperparameters
r0[i] ~ dnorm(mu.r0,tau.r0)#I(-10,10)
r1[i] ~ dnorm(mu.r1,tau.r1)#I(-10,10)
r2[i] ~ dnorm(mu.r2,tau.r2)#I(-10,10)
r3[i] ~ dnorm(mu.r3,tau.r3)#I(-10,10)
r4[i] ~ dnorm(mu.r4,tau.r4)#I(-10,10)
r5[i] ~ dnorm(mu.r5,tau.r5)#I(-10,10)
r6[i] ~ dnorm(mu.r6,tau.r6)#I(-10,10)
r7[i] ~ dnorm(mu.r7,tau.r7)#I(-10,10)
# indicator variable whether each species is exposed to sampling or not
w[i] ~ dbern(psi)
#likelihood - Ecological model for latent abundance of species i in sites j
for(j in 1:nSites){
# population abundances.
log(lambda[j,i]) <- a0[i] + a1[i]*Hyd[j] + a2[i]*Flo[j] + a3[i]*Def[j] + a4[i]*UC[j] + a5[i]*HP.C[j] + a6[i]*HP.L[j] + acom[Comunidade[j],i]
Z[j,i] ~ dpois(lambda[j,i]) # latent abundance of each species in each site
A[j,i] <- Z[j,i] * w[i] # latent abundance only for extant species
o[j,i] <- step(A[j,i]-1) # occupancy of each species in each site
# detection process model
r[j,i] <- 1/(1+exp(-(r0[i] + r1[i]*Hyd[j] + r2[i]*Flo[j] + r3[i]*Def[j] + r4[i]*UC[j] + r5[i]*HP.C[j] + r6[i]*HP.L[j] + r7[i]*Eff.2[j])))
#logit(r[j,i]) <- -(r0[i] + r1[i]*Hyd[j] + r2[i]*Flo[j] + r3[i]*Def[j] + r4[i]*UC[j] + r5[i]*HP.C[j] + r6[i]*HP.L[j] + r7[i]*Eff.2[j])
p[j,i] <- 1-pow(1-r[j,i],A[j,i])
y[j,i] ~ dbin(p[j,i], k[j]) # model observation data as binomial outcome with prob p and k trials
}#j
}#i
## counting species richness at site
#for(j in 1:nSites){
#SR[j] <- sum(o[j,]) # whole species
## counting abundance at site
#AB[j] <- sum(A[j,]) # whole species
#}
### counting total number of species
n0 <- sum(w[(nspecies+1):(nspecies+nzeros)])
N <- nspecies + n0
# Determine point level richness
for(j in 1:nSites){
Nsite[j] <- inprod(o[j,1:(nspecies+nzeros)],w[1:(nspecies+nzeros)])
}
}
",fill=TRUE)
sink()
我的数据包含 720 个站点和 29 个检测到的物种和 10 个任意未检测到的物种(用于数据增强)。
我的缩写是:
# initial values
inits < -function() list (
mu.a0=rnorm(1), mu.a1=rnorm(1), mu.a2=rnorm(1), mu.a3=rnorm(1), mu.a4=rnorm(1), mu.a5=rnorm(1), mu.a6=rnorm(1), mu.com=rnorm(1),
mu.r0=rnorm(1), mu.r1=rnorm(1), mu.r2=rnorm(1), mu.r3=rnorm(1), mu.r4=rnorm(1), mu.r5=rnorm(1), mu.r6=rnorm(1), mu.r7=rnorm(1),
sigma.a0=runif(n=1,min=0.5,max=1), sigma.a1=runif(n=1,min=0.5,max=1), sigma.a2=runif(n=1,min=0.5,max=1), sigma.a3=runif(n=1,min=0.5,max=1), sigma.a4=runif(n=1,min=0.5,max=1), sigma.a5=runif(n=1,min=0.5,max=1), sigma.a6=runif(n=1,min=0.5,max=1), sigma.com=runif(n=1,min=0.5,max=1),
sigma.r0=runif(n=1,min=0.5,max=1), sigma.r1=runif(n=1,min=0.5,max=1), sigma.r2=runif(n=1,min=0.5,max=1), sigma.r3=runif(n=1,min=0.5,max=1), sigma.r4=runif(n=1,min=0.5,max=1), sigma.r5=runif(n=1,min=0.5,max=1), sigma.r6=runif(n=1,min=0.5,max=1), sigma.r7=runif(n=1,min=0.5,max=1),
Z=matrix(rpois(n=(nspecies+nzeros)*nSites,lambda=0.05), nrow=nSites,ncol=(nspecies+nzeros)),
#Z=matrix(rbinom((nspecies+nzeros)*nSites, size=1,prob=0.1), nrow=nSites, ncol=(nspecies+nzeros)),
psi=runif(1),w=rbinom((nspecies+nzeros),1,0.5))