Jolly-Seber與WinBUGS一起運行，不與Jags

Question

我可以使用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) 

Answer 1

我想我已經想出了一個解決方案。 Jags對初始值很敏感。人們用Jags來解決這個問題的策略之一是嘗試提出接近真實值的初始值。

我用z矩陣做了這個，通過複製檢測矩陣併爲每個人從第一個檢測到最後一個檢測以1填充該拷貝。當我使用新矩陣作爲z的初始值時，運行了Jags。

以下是使用佔位方法使用WinBUGS和Jags分析樣本數據集的整個代碼。下面進一步顯示了多態方法和試驗性超羣（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運行模型時我用於爲z矩陣用於多狀態喬利-Seber模型創建初始值的代碼。我沒有再現了多態模型代碼：

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) 

到目前爲止，我已經能夠獲得超級羣（POPAN）的方法來在Jags運行的唯一方法是在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) )}