2017-07-19 4 views
2

WinBUGSで単純なJolly-Seberモデルを実行できますが、Jagsでは実行できません。私はJagsで線形回帰を実行することができ、RJagsを見つけて実行することができます。したがって、私は、Jagsがモデルコード内の1つ(またはそれ以上)の行を解釈できない可能性があると考えています。以下のコードを調べて、Jagsで実行するように修正する方法をお勧めします。当初私はおそらくprod関数がJagsで利用できないと思われました。しかし、Jagsマニュアルを検索すると、Jagsにはprod機能が含まれています。Jolly-Seberはジャグ付きではなくWinBUGSで動作します

これは私が考えることができる最も簡単な実例ですが、可能であればさらに簡素化します。

例のデータセットが下部にあります。モデルコードは、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) 
+0

通常、ジャグのエラーメッセージはかなり有益です。これを実行しようとしたときに何を伝えましたか? –

答えて

2

私は解決策が出ていると思います。 Jagsは初期値に敏感です。 Jagsでこれを克服するために戦術的な人が取るのは、真の値に近似する開始値を考え出すことです。

私はマトリックスでこれを行いました。検出マトリックスをコピーし、そのコピーを各個体の最初の検出から最後の検出までの1で満たすことによってこれを行いました。 zの初期値としてその新しいマトリックスを使用したときにJagsが実行されました。

ここでは、WinBUGSJagsの占有手法を使用してサンプルデータセットを分析するコード全体を示します。複数州アプローチと仮のスーパー人口(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) 

IはJagsでモデルを実行するときに多状態ジョリー-Seberモデルのz行列の初期値を作成するために使用されるコードです。 Iマルチステートモデルコードを再現していない:

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で実行するスーパー人口(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) )} 
関連する問題