r - Jolly-Seber 与 WinBUGS 一起运行，而不是与 Jags 一起运行

标签 r jags winbugs

我可以使用 WinBUGS 运行简单的 Jolly-Seber 模型，但不能使用 Jags 运行。我可以使用 Jags 运行线性回归，这表明 R 能够定位并执行 Jags。因此，我怀疑问题可能是 Jags 无法解释模型代码中的一行(或多行)。请检查下面的代码并建议如何修改它以在 Jags 中运行。最初我怀疑 prod 功能可能在 Jags 中不可用。但是，搜索 Jags 手册表明 Jags 确实包含 prod 函数。

这是我能想到的最简单的工作示例，但如果可能的话会进一步简化。

底部提供了示例数据集。模型代码修改自 Kery 和 Schaub (2012)。

# BUGS code

sink("C:/Users/mmiller/Documents/simple R programs/my.model.txt")
cat("
model {
for (i in 1:M) {
   for (t in 1:(n.occasions-1)) {
      phi[i,t] <- mean.phi
   }
   for (t in 1:n.occasions) {
      p[i,t] <- mean.p
   }
}
mean.phi ~ dunif(0, 1)
mean.p   ~ dunif(0, 1)
for (t in 1:n.occasions) {
   gamma[t] ~ dunif(0, 1)
}
for (i in 1:M) {
   z[i,1] ~ dbern(gamma[1])
   mu1[i] <- z[i,1] * p[i,1]
   y[i,1] ~ dbern(mu1[i])

   for (t in 2:n.occasions) {
      q[i,t-1] <- 1-z[i,t-1]
      mu2[i,t] <- phi[i,t-1] * z[i,t-1] + gamma[t] * prod(q[i,1:(t-1)])
      z[i,t] ~ dbern(mu2[i,t])
      mu3[i,t] <- z[i,t] * p[i,t]
      y[i,t] ~ dbern(mu3[i,t])
   }
}
}
",fill=TRUE)
sink()

# run R2WinBUGS

setwd('C:/Users/mmiller/Documents/simple R programs')

library(R2WinBUGS)
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])
inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.data)}  
parameters <- c("mean.p", "mean.phi")
bugs.out <- bugs(data, inits, parameters, 
                 "C:/Users/mmiller/Documents/simple R programs/my.model.txt", 
                 n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, debug=FALSE, working.dir=getwd())

print(bugs.out, digits=2)

# run R2jags

library('R2jags')
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])
inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.data)}  
parameters <- c("mean.p", "mean.phi")

jags.out2 <- jags(data = data, inits = inits, parameters, 
                  model.file = "C:/Users/mmiller/Documents/simple R programs/my.model.txt",
                  n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, working.dir=getwd())
print(jags.out2, digits=2)

这里是示例数据集:

my.data <- read.table(text = '

    1    1    1    0    0    0    0
    0    1    1    1    1    0    0
    0    1    0    0    0    0    0
    0    0    1    1    0    1    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    1    1    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    0    0    0    0    1    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    1    0    1    0    0    0    0
    0    1    0    0    0    0    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    0    0
    0    1    0    0    0    0    0
    0    0    0    0    1    1    0
    0    1    1    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    1    0    0    0
    0    0    0    1    0    0    0
    0    0    1    1    0    0    0
    0    0    1    1    1    0    0
    0    0    1    1    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    1    1
    0    0    0    0    1    1    0
    0    0    0    1    0    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    1    1
    0    0    0    1    0    0    0
    0    0    0    1    0    1    1
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    1    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1', header = FALSE)

nz <- 300
my.data <- rbind(my.data, matrix(0, ncol = ncol(my.data), nrow = nz))
dim(my.data)
head(my.data)
my.data <- as.matrix(my.data)

这是在 Jags 中运行的线性回归:

# Linear regression in JAGS using R2jags

library('R2jags')

x <- rnorm(10)
mu <- -3.2 + 1.5 * x
y <- rnorm(10, mu, sd = 4)

cat ("
model {
  for (i in 1:10) {
     y[i] ~ dnorm(mu[i], tau)
     mu[i] <- beta0 + beta1*x[i]
   }

  beta0 ~ dnorm(0, .01)
  beta1 ~ dnorm(0, .01)
  sigma ~ dunif(0,100)
  tau <- 1 / (sigma * sigma)

}

", file = 'C:/Users/mmiller/Documents/simple R programs/normal.txt')

data <- list(y=y, x=x)

inits <- function()

list(beta1 = rnorm(1), beta0 = rnorm(1), sigma = runif(1,0,2))

parameters <- c("beta0", "beta1", "sigma", "tau")

out <- jags(data = data, inits = inits, parameters, 
            model.file = 'C:/Users/mmiller/Documents/simple R programs/normal.txt',
            n.thin=1, n.chains=2, n.burnin=2000, n.iter=6000, working.dir=getwd())

print(out, digits=2)

最佳答案

我想我已经找到解决方案了。 Jags 对初始值敏感。人们使用 Jags 克服这个问题的一种策略是尝试提出接近真实值的起始值。

我使用 z 矩阵来完成此操作，方法是复制检测矩阵，并在每个个体的第一次检测到最后一次检测期间用 1 填充该副本。当我使用新矩阵作为 z 的初始值时，Jags 运行。

以下是使用 WinBUGS 和 Jags 占用方法分析示例数据集的完整代码。下面进一步展示了多状态方法和暂定 super 种群 (POPAN) 方法。

# data

my.data <- read.table(text = '

    1    1    1    0    0    0    0
    0    1    1    1    1    0    0
    0    1    0    0    0    0    0
    0    0    1    1    0    1    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    0    0    0    0
    1    0    0    0    0    0    0
    1    0    1    1    1    0    0
    1    0    0    0    0    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    0    0    0    0    1    0    0
    1    0    0    0    0    0    0
    1    1    0    0    0    0    0
    1    0    1    0    0    0    0
    0    1    0    0    0    0    0
    0    1    1    0    0    0    0
    1    0    0    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    0    0
    0    1    0    0    0    0    0
    0    0    0    0    1    1    0
    0    1    1    0    0    0    0
    0    1    1    0    0    0    0
    0    0    1    1    0    0    0
    0    0    0    1    0    0    0
    0    0    1    1    0    0    0
    0    0    1    1    1    0    0
    0    0    1    1    0    0    0
    0    0    1    0    0    0    0
    0    0    0    1    0    1    1
    0    0    0    0    1    1    0
    0    0    0    1    0    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    1    1
    0    0    0    1    0    0    0
    0    0    0    1    0    1    1
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    1    0    0    0
    0    0    0    0    1    0    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    1    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    1    0
    0    0    0    0    0    1    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    1    1
    0    0    0    0    0    1    0
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1
    0    0    0    0    0    0    1', header = FALSE)




#####
#####


my.z.init <- my.data


my.Sums <- rowSums(my.z.init)
mean(my.Sums)


first.one <- apply(my.z.init[,1:ncol(my.data)], 1, function(x) min(which(x == 1)))
first.one

last.one  <- apply(my.z.init[,1:ncol(my.data)], 1, function(x) max(which(x == 1)))
last.one


for(i in 1:nrow(my.z.init)) {

     my.z.init[i, c(first.one[i]:last.one[i])] = 1

}


my.z.init


#####
#####


nz <- 300
my.data <- rbind(my.data, matrix(0, ncol = ncol(my.data), nrow = nz))
dim(my.data)
head(my.data)
my.data <- as.matrix(my.data)

my.z.init <- rbind(my.z.init, matrix(0, ncol = ncol(my.z.init), nrow = nz))
dim(my.z.init)
head(my.z.init)
my.z.init <- as.matrix(my.z.init)


#####
#####



# BUGS code

sink("C:/Users/mmiller/Documents/simple R programs/my.model.txt")
cat("
model {
for (i in 1:M) {
   for (t in 1:(n.occasions-1)) {
      phi[i,t] <- mean.phi
   }
   for (t in 1:n.occasions) {
      p[i,t] <- mean.p
   }
}
mean.phi ~ dunif(0, 1)
mean.p   ~ dunif(0, 1)
for (t in 1:n.occasions) {
   gamma[t] ~ dunif(0, 1)
}
for (i in 1:M) {
   z[i,1] ~ dbern(gamma[1])
   mu1[i] <- z[i,1] * p[i,1]
   y[i,1] ~ dbern(mu1[i])

   for (t in 2:n.occasions) {
      q[i,t-1] <- 1-z[i,t-1]
      mu2[i,t] <- phi[i,t-1] * z[i,t-1] + gamma[t] * prod(q[i,1:(t-1)])
      z[i,t] ~ dbern(mu2[i,t])
      mu3[i,t] <- z[i,t] * p[i,t]
      y[i,t] ~ dbern(mu3[i,t])
   }
}
}
",fill=TRUE)
sink()

# run R2WinBUGS

setwd('C:/Users/mmiller/Documents/simple R programs')

library(R2WinBUGS)
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])
inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.data)}  
parameters <- c("mean.p", "mean.phi")
bugs.out <- bugs(data, inits, parameters, 
                 "C:/Users/mmiller/Documents/simple R programs/my.model.txt", 
                 n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, debug=FALSE, working.dir=getwd())

print(bugs.out, digits=2)


# run R2jags

library('R2jags')
data <- list(y = my.data, n.occasions = dim(my.data)[2], M = dim(my.data)[1])

inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), z = my.z.init)}

parameters <- c("mean.p", "mean.phi", "gamma")

jags.out2 <- jags(data = data, inits = inits, parameters, 
                  model.file = "C:/Users/mmiller/Documents/simple R programs/my.model.txt",
                  n.thin=1, n.chains=2, n.burnin=500, n.iter=1000, working.dir=getwd())
print(jags.out2, digits=2)

以下是我在 Jags 中运行模型时用于为多状态 Jolly-Seber 模型创建 z 矩阵初始值的代码。我没有重现多状态模型代码:

CH.du <- cbind(rep(0, dim(CH)[1]), CH)

my.z.init <- CH.du

first.one <- apply(my.z.init[,1:ncol(CH.du)], 1, function(x) min(which(x == 1)))
last.one  <- apply(my.z.init[,1:ncol(CH.du)], 1, function(x) max(which(x == 1)))

for(i in 1:nrow(my.z.init)) {
                                        my.z.init[i,     first.one[i]  : last.one[i]        ] = 2
     if(first.one[i] > 1)               my.z.init[i,                1  : (first.one[i] - 1) ] = 1
     if(last.one[i]  < ncol(my.z.init)) my.z.init[i, (last.one[i] + 1) : ncol(my.z.init)    ] = 3
}

nz <- 500

CH.ms <- rbind(CH.du, matrix(0, ncol = dim(CH.du)[2], nrow = nz))

CH.ms[CH.ms==0] <- 2

my.z.init.ms <- rbind(my.z.init, matrix(0, ncol = dim(my.z.init)[2], nrow = nz))

my.z.init.ms[my.z.init.ms==0] <- 1


library('jagsUI')

data <- list(y = CH.ms, n.occasions = dim(CH.ms)[2], M = dim(CH.ms)[1])

inits <- function() {list(mean.phi = runif(1, 0, 1), mean.p = runif(1, 0, 1), 

              z = cbind(rep(NA, dim(my.z.init.ms)[1]), my.z.init.ms[,-1]))}    

parameters <- c("mean.p", "mean.phi", "b", "Nsuper", "N", "B")

ni <- 2000
nt <- 3
nb <- 500
nc <- 3

js.occ <- jags(data = data, inits = inits, parameters, 
               model.file = "C:/Users/mmiller/Documents/simple R programs/js-ms.txt", 
               n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)

print(js.occ, digits = 3)

到目前为止，我能够在 Jags 中运行 super 种群 (POPAN) 方法的唯一方法是首先在 WinBUGS 中运行模型，然后使用这些平均估计值作为 Jags 的起始值。我不确定这是否是创建初始值的好方法。

inits <- function() {list(mean.phi = js.super$mean$mean.phi, 
                          mean.p   = js.super$mean$mean.p, 
                          psi      = js.super$mean$psi,
                          z        = round(js.super$mean$z)  )}

关于r - Jolly-Seber 与 WinBUGS 一起运行，而不是与 Jags 一起运行，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/45202871/

上一篇：scala - 如何将一组 RelationalGroupedDataset 传递给函数？

下一篇：django - 你能在 django 中引用多个基本模板吗？

相关文章：

R - 扩展多维数组

python - 计算R中圆半径的倒数

r - 使用 JAGS 的贝叶斯广义泊松零一技巧

r - dinterval() 用于区间删失数据？

winbugs14 - 解释 WinBUGS 陷阱以及如何使程序自动化？

r - 如何根据定义的组为树状图的标签着色？ (在 R 中)

r - 如何使用 R 中的 kableExtra::kable() 函数在标题中正确显示脚注标记

r - 将参数向量乘以 JAGS 中的自变量矩阵

r - 通过 R 在 WinBUGS 中设置种子

r - OpenBUGS 无法收敛于 WinBUGS 中收敛的模型。精度极限？

©2024 IT工具网  联系我们