openBUGS 모델을 JAGS로 변환하는 중 오류가 발생했습니다

Question

cox regression model에 대한 베타를 OpenBUGS code으로 추정하고 싶습니다. 예를 들어 베타 버전이 하나 뿐인 예제에서 수정 된 예제 코드는 다양한 베타 버전이 필요합니다. 모델에 따라 다릅니다.

이것은 내 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

두 가지가 필요합니다 모델 코드 수정 :

1) 상단의 데이터 변환은 별도의 데이터로 수행해야합니다 {}는 JAGS에서 차단됩니다 (이것은 노드 dN의 재정의에 대한 오류를 제공함).

2) 루프 인덱스 :

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

하지만 그것이 있어야 같다 :

for(k in 2:p+1){ 

는 (기인) 연산자 우선 순위에 동일하다

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

함께 이 두 가지 수정 사항은 다음 모델 코드가 테스트 코드와 함께 작동합니다.

도움이 1,363,210

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) 
} 
}" 

희망,

매트