「リスト」パッケージを使用して、ベイジアン フレームワークでリスト長分析 (基本的に多重ロジスティック回帰) を実行しようとしています。これは、ここからのみインストールできるベータ パッケージです: 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