我正在尝试在 JAGS 中建立一个允许时变协变量的生存模型。我希望它是一个参数模型——例如，假设生存遵循 Weibull 分布（但我想允许风险变化，所以指数太简单了）。因此，这本质上是可以在flexsurv包中完成的贝叶斯版本，它允许参数模型中的时变协变量。
因此，我希望能够以“计数过程”形式输入数据，其中每个主题都有多行，每行对应于其协变量保持不变的时间间隔（如本 pdf或此处所述。这是或包装允许的(start, stop]配方。survivalflexurv
不幸的是，关于如何在 JAGS 中进行生存分析的每一种解释似乎都假设每个主题只有一行。
我尝试采用这种更简单的方法并将其扩展到计数过程格式，但该模型无法正确估计分布。
失败的尝试：
这是一个例子。首先我们生成一些数据：
library('dplyr')
library('survival')
## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2
true_shape <- 2
true_scale <- 365
dat <- data_frame(person = 1:n_sub,
true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
person_start_time = runif(n_sub, min= 0, max= true_scale*2),
person_censored = (person_start_time + true_duration) > current_date,
person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)
person person_start_time person_censored person_duration
(int) (dbl) (lgl) (dbl)
1 1 11.81416 FALSE 487.4553
2 2 114.20900 FALSE 168.7674
3 3 75.34220 FALSE 356.6298
4 4 339.98225 FALSE 385.5119
5 5 389.23357 FALSE 259.9791
6 6 253.71067 FALSE 259.0032
7 7 419.52305 TRUE 310.4770
然后我们将数据分成每个受试者 2 个观察值。我只是在时间 = 300 时拆分每个主题（除非他们没有达到时间 = 300，他们只得到一个观察结果）。
## Split into multiple observations per person: --------
cens_point <- 300 # <----- try changing to 0 for no split; if so, model correctly estimates
dat_split <- dat %>%
group_by(person) %>%
do(data_frame(
split = ifelse(.$person_duration > cens_point, cens_point, .$person_duration),
START = c(0, split[1]),
END = c(split[1], .$person_duration),
TINTERVAL = c(split[1], .$person_duration - split[1]),
CENS = c(ifelse(.$person_duration > cens_point, 1, .$person_censored), .$person_censored), # <— edited original post here due to bug; but problem still present when fixing bug
TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),
END_CENS = ifelse(CENS, NA, END)
)) %>%
filter(TINTERVAL != 0)
person split START END TINTERVAL CENS TINTERVAL_CENS
(int) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl)
1 1 300.0000 0 300.0000 300.00000 1 NA
2 1 300.0000 300 487.4553 187.45530 0 187.45530
3 2 168.7674 0 168.7674 168.76738 1 NA
4 3 300.0000 0 300.0000 300.00000 1 NA
5 3 300.0000 300 356.6298 56.62979 0 56.62979
6 4 300.0000 0 300.0000 300.00000 1 NA
现在我们可以设置 JAGS 模型了。
## Set-Up JAGS Model -------
dat_jags <- as.list(dat_split)
dat_jags$N <- length(dat_jags$TINTERVAL)
inits <- replicate(n = 2, simplify = FALSE, expr = {
list(TINTERVAL_CENS = with(dat_jags, ifelse(CENS, TINTERVAL + 1, NA)),
END_CENS = with(dat_jags, ifelse(CENS, END + 1, NA)) )
})
model_string <-
"
model {
# set priors on reparameterized version, as suggested
# here: https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/d5249e71/?limit=25#8c3b
log_a ~ dnorm(0, .001)
log(a) <- log_a
log_b ~ dnorm(0, .001)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# Estimate Subject-Durations:
CENS[i] ~ dinterval(TINTERVAL_CENS[i], TINTERVAL[i])
TINTERVAL_CENS[i] ~ dweibull( nu, lambda )
}
}
"
library('runjags')
param_monitors <- c('a', 'b', 'nu', 'lambda')
fit_jags <- run.jags(model = model_string,
burnin = 1000, sample = 1000,
monitor = param_monitors,
n.chains = 2, data = dat_jags, inits = inits)
# estimates:
fit_jags
# actual:
c(a=true_shape, b=true_scale)
根据分割点的位置，模型估计基础分布的非常不同的参数。如果数据没有被拆分到计数过程表格中，它只会获取正确的参数。似乎这不是为此类问题格式化数据的方法。
如果我遗漏了一个假设，并且我的问题与 JAGS 的关系不大，而与我如何制定问题的方式更相关，那么非常欢迎提出建议。我可能对时变协变量不能用于参数生存模型（并且只能用于 Cox 模型之类的模型感到绝望，该模型假设风险恒定且实际上并不估计潜在分布）——然而，作为我在上面提到过，flexsurvregR 中的包确实适应了(start, stop]参数模型中的公式。
如果有人知道如何用另一种语言（例如 STAN 而不是 JAGS）构建这样的模型，那也将不胜感激。
编辑：
Chris Jackson 通过电子邮件提供了一些有用的建议：
我认为这里需要 JAGS 中用于截断的 T() 构造。基本上对于每个人活着但协变量不变的时期 (t[i], t[i+1])，生存时间在时期开始时被左截断，并且可能在结尾。所以你会写类似的东西
y[i] ~ dweib(shape, scale[i])T(t[i], )
我尝试按如下方式实施此建议：
model {
# same as before
log_a ~ dnorm(0, .01)
log(a) <- log_a
log_b ~ dnorm(0, .01)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# modified to include left-truncation
CENS[i] ~ dinterval(END_CENS[i], END[i])
END_CENS[i] ~ dweibull( nu, lambda )T(START[i],)
}
}
不幸的是，这并不能解决问题。使用旧代码，模型主要是正确地获得了比例参数，但在形状参数上做得非常糟糕。使用这个新代码，它非常接近正确的形状参数，但始终高估了比例参数。我注意到，高估的程度与分裂点出现的时间有关。如果分割点早（cens_point = 50），则实际上没有任何高估；如果来晚了（cens_point = 350），那么有很多。
我想这个问题可能与“重复计算”观察结果有关：如果我们在 t=300 看到一个删失的观察结果，那么来自同一个人，在 t=400 的一个未经删失的观察，对我来说，这个人似乎很直观为我们对 Weibull 参数的推断贡献了两个数据点，而实际上它们应该只贡献一个点。因此，我尝试为每个人加入随机效应；然而，这完全失败了，参数估计值很大（在 50-90 范围内）nu。我不确定为什么会这样，但也许这是一个单独的帖子的问题。由于我不是问题是否相关，您可以在此处找到整篇文章的代码，包括该模型的 JAGS代码。