Давно я пытался найти способ в rjags написать код для модели анализа выживания Байесовского Вейбулла с изменяющимся во времени зависящие от времени) ковариаты .
Данные, над которыми я работаю, касаются продолжительности от покупки до утилизации. Некоторые записи подвергаются цензуре справа. Кроме того, у меня есть несколько поясняющих переменных для каждой записи, например age
, marital status
, education level
и income
.
У меня есть код, который не содержит зависящих от времени ковариатов. Этот код работает правильно, и его результаты разумны. Теперь я хочу улучшить его, включив в него изменяющиеся во времени ковариаты, такие как age
. Вот мой простой код.
mod_string_wei <- "model{
for(i in 1 : nn) {
is.censored[i] ~ dinterval(tt[i], ee[i])
tt[i] ~ dweib(alpha,lamda[i])
lamda[i] <- exp(beta[1] + beta[2]*cc[i])
}
for (i in 1:2) {
beta[i] ~ dnorm(0.0, 1.0/1.0e5)
}
alpha ~ dgamma(1,0.0001)
# Median survival time;
median0 <- log(2)/exp(beta[1])
median1 <- log(2)/exp(beta[1]+beta[2])
#mean survival time;
mean0 <- 1/exp(beta[1])
mean1 <- 1/exp(beta[1]+beta[2])
}"
#data preparation for the model
nn <- nrow(data.s)
tt <- data.s$duration
ee <- data.s$event
cc<- data.s$ageatstart
is.censored <- data.s$is.censored
#initialisation
set.seed(72)
#Parameters
params = c("beta", "alpha")
#initial values of the model
beta.mu = -8.0
beta.sigma = 5.0
beta.num = 2
alpha.shape = 3.0
alpha.rate = 4.0
inits1 = function() {
inits = list(beta = rnorm(beta.num,beta.mu,beta.sigma), alpha = rgamma(1,alpha.shape,alpha.rate))
}
############################### 3-model Run ###############################
chain.num = 3
adapt.num = 0
#initial model
mod.syd = jags.model(textConnection(mod_string_wei), data=data, inits=inits1, n.chains=chain.num, n.adapt = adapt.num)
#burn in count
burn.count = 800
update(mod.syd, burn.count)
#sample model
thin.num = 1
iteration.num = 1000
mod.syd.sim = coda.samples(model=mod.syd,
variable.names=params, n.iter=iteration.num, thin = thin.num)
############################### 4-results ###############################
#Binding results
mod.syd.csim = do.call(rbind, mod.syd.sim)
Я был бы очень признателен, если бы кто-нибудь помог мне добавить изменяющиеся во времени (зависящие от времени) ковариаты или найти для этого учебное пособие / образец кода.