将 openBUGS 模型转换为 JAGS 时出错
Getting error while translating openBUGS model into JAGS
我想估算 cox regression model, using the OpenBUGS code 的贝塔值。
我修改了示例代码,因为在示例中它只有一个 beta,我需要各种数量的 beta,具体取决于我为其提供的模型。
这是我的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)
}
}
然而,我在JAGS中将它的语法修改为运行,它给我重新定义的错误:
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)
您需要对模型代码进行两处修改:
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)
}
}"
希望对您有所帮助,
马特
我想估算 cox regression model, using the OpenBUGS code 的贝塔值。 我修改了示例代码,因为在示例中它只有一个 beta,我需要各种数量的 beta,具体取决于我为其提供的模型。
这是我的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)
}
}
然而,我在JAGS中将它的语法修改为运行,它给我重新定义的错误:
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)
您需要对模型代码进行两处修改:
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)
}
}"
希望对您有所帮助,
马特