我正在尝试使用“Liszt”包在贝叶斯框架中执行列表长度分析(基本上是多元逻辑回归)。这是一个只能从这里安装的 beta 包: http ://www.possinghamlab.org/images/LLA/Liszt_0.8-5_1.tar.gz
当我查看包时,rjags 的模型如下:
A3ModelFun <-
function(Species, Year, ListLength = rowSums(Species), sp.names = colnames(Species),
midYear = round(median(unique(Year))),
n.iter=50000, n.burnin=20000,
resultFolder = NULL) {
.folder <-
function(tag = "ListLengthAnalysis") {
paste(gsub("[^0-9]", "_", as.character(strptime(date(), "%a %b %d %H:%M:%S %Y"))), tag, sep="_")
}
#checks
#Will not work if all the years aren't supplied
if(length(Year) != nrow(Species)) {
stop("Number of years supplied not equal to number of lists") }
#When there is only one species the analysis wont run ...
if(length(Species) == 1) {
stop("Only one species") }
#Assign default folder name if unspecified
if(is.null(resultFolder)) resultFolder <- .folder()
is.folder <- function(file) file.info(file)$isdir
if(is.na(is.folder(resultFolder))) dir.create(resultFolder) else print("folder exists")
## Set up a folder for the results
# if(file.exists(resultFolder))
# if(!is.folder(resultFolder)) stop("File ", resultFolder, " exists and is not a folder") else
# dir.create(resultFolder)
#Assign default midYear if unspecified
if(is.null(midYear)) midYear <- round(median(unique(Year)))
## Create a function to generate inits (outside the for loop as it is the same params for every model)
inits <- function() list(a1=rnorm(1), a2=rnorm(1), a3=rnorm(1))
ModelFile <- tempfile()
LL <- ListLength
no <- nrow(Species) ## need to set these
my <- midYear
nyear <- length(unique(Year))
meanLogLL <- mean(log(ListLength))
halfyr <- nyear/2
##Progress Report
j <- 1
tot <- length(sp.names)
for (i in sp.names) {
cat(i, "\n")
cat(j, "out of", tot, "\n")
j <- j+1
}
#####List Length Analysis#####
for (i in sp.names) { # changed this from for (i in sp.names[]) { & for (i in sp.names[i])
## Define data for each model (we can reuse this object for each model)
LLAData.1sp <- list(pres=Species[, i], yr=Year, LL=ListLength)
## Define the JAGS model
cat(sprintf('
model {
for (i in 1:%s) {
pres[i] ~ dbern(p[i])
logit(p[i]) <- a1 +
a2 * (log(LL[i]) - %s) +
a3 * (yr[i] - %s)
}
a1 ~ dnorm(0, 0.0001)
a2 ~ dnorm(0, 0.0001)
a3 ~ dnorm(0, 0.0001)
}', no, meanLogLL, my, my),
file=ModelFile)
## Run a model and store it in an R object called model_spname (where spname changes with each model)
assign(sprintf('model_%s', i),
jags(LLAData.1sp, inits,
parameters.to.save=c(
'a1', 'a2', 'a3'),
n.iter=n.iter, n.burnin=n.burnin, ModelFile))
## Save the current model to a file with the same names
save(list=sprintf('model_%s', i), file=file.path(resultFolder, sprintf('model_%s.RData', i)))
## Delete the current model to save space
rm(list = sprintf('model_%s', i))
}
file.remove(ModelFile)
return(resultFolder)
}
模型本身有很多不相关的东西,但我觉得最好放置整段代码。
当我用我的数据运行它时,我总是得到错误:
Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, :
Error in node pres[5]
Node inconsistent with parents
我发现有更多人收到此错误消息,但它始终与与所用分布相关的数据有关(即伽马分布的负数据,二项式的泊松数据)。我有二项式分布式数据,但仍然得到错误。我不知道为什么。
我的模型如下(有关我的数据框的 20 行子样本,请参见下面的 dput,我仅在具有完整数据集的不同节点上得到相同的错误):
library(dplyr)
df <- (dput)
sp.ll <-as.matrix(df %>% select(-Year, -Hok_uniek, -LL))
mod <- A3ModelFun(sp.ll, df$Year, df$LL)
数据(在驱动器上,否则我将超过最大字符数):
https://drive.google.com/open?id=1etj_7IrEHJbYzBk9p_8tW6B-X0HpbLU9WMWMJL-QT_I