Ошибка при переводе модели 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)){ 

С эти две модификации, следующий код модели работает для меня с вашим тестовым кодом:

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

Надежда, что помогает,

Matt