Написание модели Кокфа в JAGS с использованием R (Байесовский анализ)

Question

Я новичок в Байесовском анализе, и у меня проблемы с программированием (думаю, что это проблема). Я пытаюсь выполнить байесовский анализ выживаемости, используя MCMC с JAGS. У меня есть моя простая модель записи с функцией выживания:

fit0 <- coxph(Surv(time, event) ~ group + var1 + var2 + var3, data=mydata_simple)
summary(fit0)

Для основной переменной, которая меня интересует (группа (0,1)), я получаю HR: 0,40 (р = 0,023).

Я пытаюсь провести этот анализ с байесовским подходом с использованием неинформативных априорных значений. Для этого я пишу в JAGS, используя функцию cat.

cat("
model{

for(i in 1:N){

sVL[i] <- (var1[i]) 

sCD4 [i]<- (var2[i])

h0[i]<-lambda*alpha*pow(time[i],alpha-1)

h[i]<-h0[i]*exp(beta[1]*equals(group[i],1) +
beta[2]*var3[i] + beta[3]*sVL[i] + beta[4]*sCD4[i])

#Cumulative Hazard H[t]
cumHaz[i]<-lambda*pow(time[i],alpha)*exp(beta[1]*equals(group[i],1) +
beta[2]*var3[i] + beta[3]*sVL[i]+ beta[4]*sCD4[i])
log.Surv[i] <- -cumHaz[i]
Surv[i]<-exp( -cumHaz[i])



# Definition of the log-likelihood using zeros trick
phi[i] <- C - ((1 - event[i]) * log(h[i]) + log.Surv[i])
zeros[i] ~ dpois(phi[i])
}


for(i in 1:5){

beta[i] ~ dnorm(0,0.001)
#Parameters probabilities
prob[i]<-step(beta[i])

#Hazard Ratio probabilities
HR[i]<-exp(beta[i])
probHR[i] <- step(HR[i]-1)

}
lambda~dgamma(0.1,0.1)
sigma~dunif(0,100)
alpha<-1/sigma


}
", file= "my_weibullph.txt")

Тогда я могу запустить JAGS:

ni <- 5000
na <- 500
nt <- 5
data <- list(N=111, time = mydata_simple$time, event = mydata_simple$event,group = mydata_simple$group, var3 =mydata_simple$var3, var2 = mydata_simplevar2, var1 = mydata_simple$var1, C=50000,zeros=rep(0,111))
inits <- function(){list(beta=rnorm(5,0,0.5), sigma=runif(1,0,1), lambda=runif(1,0.01,1))}
parameters <- c("beta","alpha","lambda","prob" ,"HR", "probHR","h0")

model_ph<- jags.model(data=data, file="my_weibullph.txt",inits=inits,n.adapt=na,n.chains=3)
model_ph_result <- coda.samples(model_ph, variable.names=parameters,
                                  n.iter=ni, thin=nt)#To generate posterior samples in mcmcm.list
formatsave(model_ph_result, file="model_ph_new.RData")

Чтобы получить результаты HR:

bss <- do.call(rbind,suppression_ph_result)
sims.list <- vector("list", length(parameters))
names(sims.list) <- parameters
for(p in seq_along(parameters)){
  iik <- grep(paste("^", parameters[p], sep=""), colnames(bss))
  sims.list[[p]] <- bss[,iik]
}

# Hazard ratios
HRB_group<-sims.list[[5]][,2]
HRB_var1<-sims.list[[5]][,3]
HRB_var2<-sims.list[[5]][,4]
HRB_var3<-sims.list[[5]][,5]

#prob HR
pHRB_group<-sims.list[[6]][,1]
pHRB_var1<-sims.list[[6]][,2]
pHRB_var2<-sims.list[[6]][,3]
pHRB_var<-sims.list[[6]][,4]

Хорошо, когда я получу среднее значениеHRB_group я получаю 1,18 (вместо 0,40, как в пакете с выживанием пакета). Я не знаю, что я делаю неправильно в сценарии, чтобы получить такие разные результаты ...

Заранее спасибо. Извините, если я не могу предоставить данные, но они являются конфиденциальными, и я не знаю, чтобы смоделировать другие данные, чтобы включить их здесь.