计数过程中参数化生存模型的表示JAGS格式

我试图在JAGS中建立一个生存模型，允许时变协变量。我希望它是一个参数模型——例如，假设生存遵循威布尔分布(但我希望允许风险变化，所以指数太简单了)。因此，这本质上是在flexsurv包中可以做什么的贝叶斯版本，它允许参数化模型中的时变协变量。

因此，我希望能够以“计数过程”形式输入数据，其中每个主题有多行，每一行对应于一个时间间隔，其中协变量保持不变(如在这个PDF中或在这里所述)。这是survival或flexurv包允许的(start, stop]公式。

不幸的是，每一个关于如何在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模型这样的模型，它假设恒定的危险，并且实际上不估计潜在的分布)——然而，正如我上面提到的，R中的flexsurvreg包确实适合参数化模型中的(start, stop]公式。

如果有人知道如何用另一种语言(例如STAN而不是JAGS)构建这样的模型，那也会很感激。

编辑:

克里斯·杰克逊通过电子邮件提供了一些有用的建议:

我认为这里需要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时看到一个未经审查的观察结果，对我来说，直觉上这个人为我们关于威布尔参数的推断贡献了两个数据点，而实际上他们应该只贡献一个点。因此，我尝试为每个人加入随机效应;然而，这完全失败了，对nu参数的估计很大(在50-90范围内)。我不确定这是为什么，但也许这是一个单独的问题。因为我不知道这些问题是否相关，你可以找到这篇文章的代码，包括该模型的JAGS代码在这里。