jwd*_*ink 65 r bayesian jags survival-analysis weibull
我正在尝试在JAGS中建立一个生存模型,允许时变协变量.我希望它是一个参数模型 - 例如,假设生存遵循威布尔分布(但我想让危险变化,所以指数太简单了).因此,这基本上是可以在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],)
}
}
不幸的是,这并没有完全解决问题.使用旧代码,模型主要是使缩放参数正确,但在shape参数上做得非常糟糕.使用这个新代码,它非常接近正确的形状参数,但始终高估了scale参数.我注意到过度估计的程度与分裂点的延迟时间有关.如果分裂点是早期(cens_point = 50),那么实际上没有任何过高估计; 如果它迟到(cens_point = 350),有很多.
我想也许这个问题可能与"重复计算"观察结果有关:如果我们在t = 300时看到一个被审查的观察,那么从同一个人那里,在t = 400时进行未经审查的观察,这对我来说似乎很直观我们对Weibull参数的推断贡献了两个数据点,实际上它们应该只贡献一个点.因此,我尝试为每个人加入随机效应; 然而,这完全失败了,对nu参数的估计很大(在50-90范围内).我不确定为什么会这样,但也许这是一个单独帖子的问题.因为我不是问题是否相关,你可以找到的代码,这整个的帖子,包括该模型尖齿代码,在这里.
您可以使用rstanarmpackage，它是 STAN 的包装器。它允许使用标准 R 公式符号来描述生存模型。stan_surv函数接受“计数过程”形式的参数。包括威布尔在内的不同基本危险函数可用于拟合模型。
rstanarm-功能的生存部分stan_surv在 CRAN 上仍然不可用，因此您应该直接从mc-stan.org安装该软件包。
install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))\n请看下面的代码：
\nlibrary(dplyr)\nlibrary(survival)\nlibrary(rstanarm)\n\n## Make the Data: -----\nset.seed(3)\nn_sub <- 1000\ncurrent_date <- 365*2\n\ntrue_shape <- 2\ntrue_scale <- 365\n\ndat <- data_frame(person = 1:n_sub,\n true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),\n person_start_time = runif(n_sub, min= 0, max= true_scale*2),\n person_censored = (person_start_time + true_duration) > current_date,\n person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)\n)\n\n## Split into multiple observations per person: --------\ncens_point <- 300 # <----- try changing to 0 for no split; if so, model correctly estimates\ndat_split <- dat %>%\n group_by(person) %>%\n do(data_frame(\n split = ifelse(.$person_duration > cens_point, cens_point, .$person_duration),\n START = c(0, split[1]),\n END = c(split[1], .$person_duration),\n TINTERVAL = c(split[1], .$person_duration - split[1]), \n CENS = c(ifelse(.$person_duration > cens_point, 1, .$person_censored), .$person_censored), # <\xe2\x80\x94 edited original post here due to bug; but problem still present when fixing bug\n TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),\n END_CENS = ifelse(CENS, NA, END)\n )) %>%\n filter(TINTERVAL != 0)\ndat_split$CENS <- as.integer(!(dat_split$CENS))\n\n\n# Fit STAN survival model\nmod_tvc <- stan_surv(\n formula = Surv(START, END, CENS) ~ 1,\n data = dat_split,\n iter = 1000,\n chains = 2,\n basehaz = "weibull-aft")\n\n# Print fit coefficients\nmod_tvc$coefficients[2]\nunname(exp(mod_tvc$coefficients[1]))\n输出，与真实值一致（true_shape <- 2; true_scale <- 365）：
> mod_tvc$coefficients[2]\nweibull-shape \n 1.943157 \n> unname(exp(mod_tvc$coefficients[1]))\n[1] 360.6058\n您还可以查看 STAN 源代码，rstan::get_stanmodel(mod_tvc$stanfit)将 STAN 代码与您在 JAGS 中所做的尝试进行比较。