是否可以为 JAGS 中的循环定义随机限制?
Is possible to define a random limit for a loop in JAGS?
我正在尝试按照 Hui、Ibrahim 和 Sinha (1999) 概述的方法实施具有治愈分数的 Weibull 比例风险模型 - 一种新的具有生存分数的生存数据贝叶斯模型。但是,我不确定是否可以为 JAGS 中的循环定义随机限制。
library(R2OpenBUGS)
library(rjags)
set.seed(1234)
censored <- c(1, 1)
time_mod <- c(NA, NA)
time_cens <- c(5, 7)
tau <- 4
design_matrix <- rbind(c(1, 0, 0, 0), c(1, 0.2, 0.2, 0.04))
jfun <- function() {
for(i in 1:nobs) {
censored[i] ~ dinterval(time_mod[i], time_cens[i])
time_mod[i] <- ifelse(N[i] == 0, tau, min(Z))
for (k in 1:N[i]){
Z[k] ~ dweib(1, 1)
}
N[i] ~ dpois(fc[i])
fc[i] <- exp(inprod(design_matrix[i, ], beta))
}
beta[1] ~ dnorm(0, 10)
beta[2] ~ dnorm(0, 10)
beta[3] ~ dnorm(0, 10)
beta[4] ~ dnorm(0, 10)
}
inits <- function() {
time_init <- rep(NA, length(time_mod))
time_init[which(!status)] <- time_cens[which(!status)] + 1
out <- list(beta = rnorm(4, 0, 10),
time_mod = time_init,
N = rpois(length(time_mod), 5))
return(out)
}
data_base <- list('time_mod' = time_mod, 'time_cens' = time_cens,
'censored' = censored, 'design_matrix' = design_matrix,
'tau' = tau,
'nobs' = length(time_cens[!is.na(time_cens)]))
tc1 <- textConnection("jmod", "w")
write.model(jfun, tc1)
close(tc1)
# Calling JAGS
tc2 <- textConnection(jmod)
j <- jags.model(tc2,
data = data_base,
inits = inits(),
n.chains = 1,
n.adapt = 1000)
我观察到以下错误:
Error in jags.model(tc2, data = data_base, inits = inits(), n.chains = 1, :
RUNTIME ERROR:
Compilation error on line 6.
Unknown variable N
Either supply values for this variable with the data
or define it on the left hand side of a relation.
我不完全确定,但我很确定你不能在一般的 BUGS 中声明随机数量的节点,所以它不会是特定的 JAGS 的怪癖。
不过,您可以找到解决方法。
由于 BUGS 是一种声明性语言而不是过程性语言,因此声明任意但 确定性 节点数(假设 "large enough")就足够了,然后仅将 随机 个节点与分布和观察到的数据相关联,使其余节点具有确定性。
一旦观察到 N[i]
的最大值(比方说 N.max
),就可以将其作为 参数 传递给 JAGS,然后更改你的代码:
for (k in 1:N[i]){
Z[k] ~ dweib(1, 1)
}
进入这个:
for (k in 1:N.max){
if (k <= N[i]){
Z[k] ~ dweib(1, 1)
} else {
Z[k] <- 0
}
}
我希望这能解决您的问题。所以请稍后反馈。
不用说,如果您有一些 non-zero 观察到的数据与确定性 Z[k]
相关联,那么 Jags 内部的一切都会崩溃...
我正在尝试按照 Hui、Ibrahim 和 Sinha (1999) 概述的方法实施具有治愈分数的 Weibull 比例风险模型 - 一种新的具有生存分数的生存数据贝叶斯模型。但是,我不确定是否可以为 JAGS 中的循环定义随机限制。
library(R2OpenBUGS)
library(rjags)
set.seed(1234)
censored <- c(1, 1)
time_mod <- c(NA, NA)
time_cens <- c(5, 7)
tau <- 4
design_matrix <- rbind(c(1, 0, 0, 0), c(1, 0.2, 0.2, 0.04))
jfun <- function() {
for(i in 1:nobs) {
censored[i] ~ dinterval(time_mod[i], time_cens[i])
time_mod[i] <- ifelse(N[i] == 0, tau, min(Z))
for (k in 1:N[i]){
Z[k] ~ dweib(1, 1)
}
N[i] ~ dpois(fc[i])
fc[i] <- exp(inprod(design_matrix[i, ], beta))
}
beta[1] ~ dnorm(0, 10)
beta[2] ~ dnorm(0, 10)
beta[3] ~ dnorm(0, 10)
beta[4] ~ dnorm(0, 10)
}
inits <- function() {
time_init <- rep(NA, length(time_mod))
time_init[which(!status)] <- time_cens[which(!status)] + 1
out <- list(beta = rnorm(4, 0, 10),
time_mod = time_init,
N = rpois(length(time_mod), 5))
return(out)
}
data_base <- list('time_mod' = time_mod, 'time_cens' = time_cens,
'censored' = censored, 'design_matrix' = design_matrix,
'tau' = tau,
'nobs' = length(time_cens[!is.na(time_cens)]))
tc1 <- textConnection("jmod", "w")
write.model(jfun, tc1)
close(tc1)
# Calling JAGS
tc2 <- textConnection(jmod)
j <- jags.model(tc2,
data = data_base,
inits = inits(),
n.chains = 1,
n.adapt = 1000)
我观察到以下错误:
Error in jags.model(tc2, data = data_base, inits = inits(), n.chains = 1, :
RUNTIME ERROR:
Compilation error on line 6.
Unknown variable N
Either supply values for this variable with the data
or define it on the left hand side of a relation.
我不完全确定,但我很确定你不能在一般的 BUGS 中声明随机数量的节点,所以它不会是特定的 JAGS 的怪癖。
不过,您可以找到解决方法。
由于 BUGS 是一种声明性语言而不是过程性语言,因此声明任意但 确定性 节点数(假设 "large enough")就足够了,然后仅将 随机 个节点与分布和观察到的数据相关联,使其余节点具有确定性。
一旦观察到 N[i]
的最大值(比方说 N.max
),就可以将其作为 参数 传递给 JAGS,然后更改你的代码:
for (k in 1:N[i]){
Z[k] ~ dweib(1, 1)
}
进入这个:
for (k in 1:N.max){
if (k <= N[i]){
Z[k] ~ dweib(1, 1)
} else {
Z[k] <- 0
}
}
我希望这能解决您的问题。所以请稍后反馈。
不用说,如果您有一些 non-zero 观察到的数据与确定性 Z[k]
相关联,那么 Jags 内部的一切都会崩溃...