更新：我运行了 JAGS 模型，这消除了我问题中令人分心的部分。这实际上是关于为 dinterval() 和 inits 正确准备数据。
当从时间向量和删失状态向量开始时，我找不到准备右删失生存数据以用于 JAGS 的具体示例。我见过的例子都有一个恒定的审查时间。
我的两个问题是：
- 当您从事件时间向量和审查/失败状态的匹配向量开始时，准备 dinterval(t[i], c[i]) 的正确方法是什么？
- 提供事件时间的初始值向量的原因是什么，其中未审查时间设置为 NA 并且审查时间设置为大于观察审查时间的值？我真的很难理解这一点。
更新：我想我有（2）的答案。显然，如果您尝试覆盖 init 向量中的值，JAGS 会抱怨，但如果该向量被初始化为 NA 它可以工作。
请注意，我对这些要求的理解来自于查看现有的在线示例。如果要求没有根据，我很乐意学习。
我使用来自生存包的 aml 数据集创建了一个简短的测试，该数据集进行生存拟合，然后尝试使用 JAGS 拟合 Weibull 分布并添加这些点以进行比较。最初我认为您必须为 dinterval(t[i], c[i]) 的 c[i] 中的未经审查的时间传递零，但这会导致错误：
节点中的错误 is.censored[1] 在初始化时观察到的节点与观察到的父节点不一致。尝试设置适当的初始值。
这是我的示例 R 脚本。有些行似乎缩进太多，不知道为什么。我将所有选项卡转换为空格，它在我的编辑窗口中看起来很干净：
# Right-censoring test
require(survival)
require(runjags)
graphics.off()
rm(fit, S, results)
# SURVIVAL FIT AND PLOT
fit <- survfit(Surv(aml$time, aml$status) ~ 1, conf.type = "log-log")
plot(fit, xlab = "time", ylab = "survival", main = "AML survival")
# JAGS WEIBULL FIT WITH CENSORING
# Model uses three different versions of observed times.
# t.orig is the original set of observed times including censored times.
# t.mod has the censored times set to NA.
# t.cens is the set I have questions about. Currently the same as t.orig
# but is that the proper choice?
#
model <- "model {
for(i in 1:N) {
is.censored[i] ~ dinterval(t.mod[i], t.cens[i])
t.mod[i] ~ dweib(a, b)
S[i] <- 1 - pweib(t.orig[i], a, b)
}
a ~ dgamma(3.5, 2)
b ~ dgamma(1.5, 100)
}"
censored <- !aml$status
t.orig <- aml$time
t.mod <- aml$time
t.mod[censored] <- NA
t.cens <- aml$time
# Online examples seem to indicate special treatment of the c[i]
# passed to dinterval(t[i], c[i]), changing the uncensored times.
#t.cens[!censored] <- 0 # THIS CAUSES AN ERROR
data <- dump.format(list(
t.orig = t.orig,
t.mod = t.mod,
N = length(t.mod),
t.cens = t.cens,
is.censored = as.numeric(censored)
))
# Apparently you must provide initial values of the time vector
# with uncensored times of NA and censored times > the observed time.
# I don't understand this requirement at all.
t.init <- rep(NA, length(t.mod))
t.init[censored] <- t.cens[censored] + 1
inits <- c(
dump.format(list(a = 1.1, b = 3.0e-4, t.mod = t.init)),
dump.format(list(a = 2.3, b = 0.02, t.mod = t.init))
)
results = run.jags(model, monitor=c("S", "a", "b"), data=data, inits=inits)
S <- print(results, vars="S")
points(aml$time, S[,4], col="blue", pch=19)
segments(aml$time, S[,1], aml$time, S[,3], col="blue", lwd=2)
感谢您花时间看这个。