1

我正在使用 R2 Winbugs 在 Kery 等人 2009 年（此处为论文）之后使用调查的真实数据运行开放总体二项式混合模型。

当我第一次尝试时，我只使用一个协变量来表示我自己制作的丰度 (X1)。该模型编译并运行，但我得到“未定义的真实结果”。

如果可能的话，任何人都可以帮助我找到解决这个问题的方法吗？

我的计数数据可以在这里下载

代码如下

#data formatting
bdat <- read.table("PyroniaC_model_test2.txt", header = TRUE)
names(bdat)
#convert data to array
y <- array(NA, dim = c(620, 2, 352))  # 620 sites, 2 reps, 352 days
str(y)
 for(k in 1:352){
 sel.rows <- bdat$Bweekcoded == k
 y[ , , k] <- as.matrix((bdat)[sel.rows, 3:4],stringsAsFactors = F)
 }
#covs
X1<-rnorm(620, 0,100)

# Specify model in BUGS language
sink("Nmix3C.txt")
cat("
model {

# Priors
alpha0 ~ dunif(0, 0.2) #intercept log(N)
beta1 ~ dunif(0, 2) #beta cov1
r.rate ~ dunif(-1, 1) #growth rate
for (k in 1:352){   #prior for detection
 p[k] ~ dunif(0, 1)
}
# Likelihood
# Ecological model for true abundance
for (k in 1:352){                          # Loop over weeks (352)
for (i in 1:R){                       # Loop over R sites (620)
log(lambda[i,k])<-alpha0+beta1*X1[i]+r.rate*(k-1)
N[i,k] ~ dpois(lambda[i,k])          # Abundance

# Observation model for replicated counts
for (j in 1:T){                    # Loop over temporal reps (2)
y[i,j,k] ~ dbin(p[k], N[i,k])   # Detection

# Assess model fit using Chi-squared discrepancy
# Compute fit statistic E for observed data
eval[i,j,k] <- p[k] * N[i,k]   # Expected values
E[i,j,k] <- pow((y[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k] + 0.5)
# Generate replicate data and compute fit stats for them
y.new[i,j,k] ~ dbin(p[k], N[i,k])
E.new[i,j,k] <- pow((y.new[i,j,k] - eval[i,j,k]),2) / (eval[i,j,k] + 0.5)

} #j 
} #i 
} #k 

# Derived and other quantities
for (k in 1:352){
totalN[k] <- sum(N[,k])# Total pop. size across all sites
#mean.abundance[k] <- exp(alpha.lam[k])
}
fit <- sum(E[,,])
fit.new <- sum(E.new[,,])
}

",fill = TRUE)
sink()

# Bundle data
R = nrow(y)
T = ncol(y)
win.data <- list(y = y, R = R, T = T, X1=X1)

 # Initial values
Nst <- apply(y, c(1, 3), max) +1
Nst[is.na(Nst)] <- 1
inits <- function(){list(N = Nst , alpha0 = runif(1, 0, 0.2),r.rate=runif(1,-1,1),beta1 = runif(1,     0,2), p=runif(1,0,1))}

# Parameters monitored
params <- c("totalN", "alpha0", "beta1","p", "fit", "fit.new", "r.rate")

# MCMC settings
ni <- 10000
nt <- 8
nb <- 2000
nc <- 3

library(R2WinBUGS)
# Call WinBUGS from R (BRT 1 min)
out3C <- bugs(win.data, inits, params, "Nmix3C.txt", n.chains = nc, n.thin = nt, n.iter = ni,  n.burnin = nb, debug = TRUE)

非常感谢你们所有人！

PS。我也尝试了一个真正的协变量 X1，但我得到了同样的错误。

4

0 回答 0