openBUGSモデルをJAGSに翻訳する際にエラーが発生する

Question

cox regression modelのベータ版をOpenBUGS codeで見積もりたいとします。 例のコードを変更しました。例ではベータが1つしかないため、ベータ版の数はさまざまです。モデルに依存しています。

これは私のopenBUGSモデルです。期待どおりに動作します：

bugsmodel <- function(){ 
# Set up data 
for(i in 1:N) { 
    for(j in 1:T) { 
    Y[i,j] <- step(obs.t[i] - t[j] + eps) 
    dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i] 
    } 
} 
# Model 
for(i in 1:N){ 
    betax[i,1] <- 0 
    for(k in 2:p+1){ 
    betax[i,k] <- betax[i,k-1] + beta[k-1]*x[i,k-1] 
    } 
} 
for(j in 1:T) { 
    for(i in 1:N) { 
    dN[i, j] ~ dpois(Idt[i, j]) # Likelihood 
    Idt[i, j] <- Y[i, j] * exp(betax[i,p+1]) * dL0[j] # Intensity 
    } 
    dL0[j] ~ dgamma(mu[j], c) 
    mu[j] <- dL0.star[j] * C# prior mean hazard 
} 
c <- 0.001 
r <- 0.1 
for (j in 1 : T) { 
    dL0.star[j] <- r * (t[j + 1] - t[j]) 
} 
for(k in 1:p){ 
    beta[k] ~ dnorm(0.0,0.000001) 
} 
} 

はしかし、私はそれが私の再定義のエラーを与え、ぎざぎざで実行するようにその構文を修正：

model_jags <- "model{ 
    # Set up data 
for(i in 1:N) { 
    for(j in 1:T) { 
    Y[i,j] <- step(obs.t[i] - t[j] + eps) 
    dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i] 
    } 
} 
# Model 
for(i in 1:N){ 
    betax[i,1] <- 0 
    for(k in 2:p+1){ 
    betax[i,k] <- betax[i,k-1] + beta[k-1]*x[i,k-1] 
    } 
} 
for(j in 1:T) { 
    for(i in 1:N) { 
    dN[i, j] ~ dpois(Idt[i, j]) # Likelihood 
    Idt[i, j] <- Y[i, j] * exp(betax[i,p+1]) * dL0[j] # Intensity 
    } 
    dL0[j] ~ dgamma(mu[j], c) 
    mu[j] <- dL0.star[j] * C# prior mean hazard 
} 
c <- 0.001 
r <- 0.1 
for (j in 1 : T) { 
    dL0.star[j] <- r * (t[j + 1] - t[j]) 
} 
for(k in 1:p){ 
    beta[k] ~ dnorm(0.0,0.000001) 
} 
}" 

テストコード：

あなたが必要

n = 100 
round=2 
x1 = rbinom(n,size=1,prob=0.5) 
x2 = rbinom(n,size=1,prob=0.5) 
x3 = rbinom(n,size=1,prob=0.5) 
x = t(rbind(x1,x2,x3)) 
censortime = runif(n,0,1) 
survtime= rexp(n,rate=exp(x1+2*x2+3*x3)) 
survtime = round(survtime,digits=round) 
event = as.numeric(censortime>survtime) 
y = survtime; 
y[event==0] = censortime[event==0] 
t=sort(unique(y[event==1])) 
t=c(t,max(censortime)) 
bigt=length(t)-1 
##################################### 
model=c(1,1,1) 
x <- x[,model==1] 
p <- sum(model) # models have betas with 1 
params <- c("beta","dL0") 
data <- list(x=x,obs.t=y,t=t,T=bigt,N=n,fail=event,eps=1E-10,p=p) 
inits <- function(){list(beta = rep(0,p), dL0 = rep(0.0001,bigt))} 

jags <- jags.model(textConnection(model_jags), 
       data = data, 
       n.chains = 1, 
       n.adapt = 100) 

Answer 1

2モデルコードの変更：

1）先頭のデータ変換は別のデータで行う必要があります{}はJAGSでブロックします（これはノードdNの再定義に関するエラーを出しています）。

2）ループインデックス：

for(k in (2:p)+1){ 

しかし、私はそれがあるべき推測：

for(k in 2:p+1){ 

は（原因）演算子の優先順位と同じです

for(k in 2:(p+1)){ 

付きこれらの2つの変更は、あなたのテストコードで私のために次のモデルコードが機能します：

助け

model_jags <- " 
data{ 
    # Set up data 
for(i in 1:N) { 
    for(j in 1:T) { 
    Y[i,j] <- step(obs.t[i] - t[j] + eps) 
    dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * fail[i] 
    } 
} 
} 

# Model 
model{ 
for(i in 1:N){ 
    betax[i,1] <- 0 
    for(k in 2:(p+1)){ 
    betax[i,k] <- betax[i,k-1] + beta[k-1]*x[i,k-1] 
    } 
} 
for(j in 1:T) { 
    for(i in 1:N) { 
    dN[i, j] ~ dpois(Idt[i, j]) # Likelihood 
    Idt[i, j] <- Y[i, j] * exp(betax[i,p+1]) * dL0[j] # Intensity 
    } 
    dL0[j] ~ dgamma(mu[j], c) 
    mu[j] <- dL0.star[j] * C# prior mean hazard 
} 
c <- 0.001 
r <- 0.1 
for (j in 1 : T) { 
    dL0.star[j] <- r * (t[j + 1] - t[j]) 
} 
for(k in 1:p){ 
    beta[k] ~ dnorm(0.0,0.000001) 
} 
}" 

希望、

マット