survival-analysis - JAGS 中的逐步右删失生存分析

Question

这是对 SE 早期帖子的一种跟进:https://stats.stackexchange.com/questions/70858/right-censored-survival-fit-with-jags

但是在这里，我希望看到一个完整的 R 脚本(从头到尾)对 JAGS 中的右删失数据进行生存分析。我发现的所有站点都需要非常熟练地使用 JAGS，因此我很难理解如何从一行代码转到另一行代码。我知道有很多问题要问......

无论如何，这里有一些示例生存数据。组是 t1、t2、t3。 NA 指的是右删失数据(删失截止值 = 3)。

t1 <- c(1.73, NA, NA, NA, NA,0.77, NA, NA, NA, NA, NA, NA,0.5,1.06, NA, NA, NA, NA, NA,0.46, NA)
t2 <- c(1.42, NA, NA, NA, NA, NA,0.69,1.84, NA, NA, NA,1.47,1.6,1.8, NA, NA, NA, NA, NA,0.73, NA,1.28,3,2.97)
t3 <- c(0.88, NA, NA,1.65,1.75, NA, NA,1.01,1.46,1.95, NA, NA,2.93, NA,0.78,1.05,1.52, NA)


#Specify model in BUGS language


sink("model.txt")
cat("
model
{




}
",fill = TRUE)
sink()

#Bundle data
data<- list()

#Parameters monitored
parameters<-c()

#Initial values
inits <- list(

# MCMC settings
ni <-  
nt <- 
nb <- 
nc <- 


fit <- jags(data, inits, parameters, "model.txt", n.iter=ni, n.thin=nt, n.burnin=nb, n.chains=nc, working.directory = getwd())

我知道有很多问题要问，但我花了几天时间试图将一些东西拼凑起来，但我一直迷路/困惑。我知道现在有一些包可以运行这种分析，但我真的很想学习如何从头开始并自己构建它!谢谢各位读者!

Answer 1

我没有做很多生存分析(并且你没有说明你想在这部分使用哪个分布 - 有不同的选项)，但是这段代码应该让你开始进行区间审查:

    library("runjags")

    # Your data
    t1 <- c(1.73, NA, NA, NA, NA,0.77, NA, NA, NA, NA, NA, NA,0.5,1.06, NA, NA, NA, NA, NA,0.46, NA)
    t2 <- c(1.42, NA, NA, NA, NA, NA,0.69,1.84, NA, NA, NA,1.47,1.6,1.8, NA, NA, NA, NA, NA,0.73, NA,1.28,3,2.97)
    t3 <- c(0.88, NA, NA,1.65,1.75, NA, NA,1.01,1.46,1.95, NA, NA,2.93, NA,0.78,1.05,1.52, NA)

    # Combine these into a single vector to make the code cleaner
    alldata <- rbind(cbind(t1, 1), cbind(t2, 2), cbind(t3, 3))
    T.obs <- alldata[,1]
    Group <- alldata[,2]
    N <- length(T.obs)

    # The censoring setup - in this case 0 for T.obs < 3 and 1 for T.obs > 3
    Censoring <- as.numeric(is.na(T.obs))
    Breakpoint <- 3

    # A simple JAGS model:
    m <- "
    model{
        for(i in 1:N){
            # The censoring part:
            Censoring[i] ~ dinterval(T.obs[i], Breakpoint)
            # The regression part - you may well want to change dexp to something else:
            T.obs[i] ~ dexp(rate[Group[i]])
        }   
        rate[1] ~ dgamma(0.01, 0.01)
        rate[2] ~ dgamma(0.01, 0.01)
        rate[3] ~ dgamma(0.01, 0.01)

        #data# N, Censoring, Breakpoint, T.obs, Group
        #monitor# rate, T.obs
    }
    "

    # One of the things we need to do is help JAGS initialise T.obs:
    T.obs.init <- ifelse(is.na(T.obs), 4, NA)

    # The function call:
    results <- run.jags(m, n.chains=2, inits=list(T.obs=T.obs.init))

    # Look at results:
    results

这使用 runjags 包，它执行一些自动收敛等诊断，并允许在模型代码而不是 R 代码中速记使用 #data# 和 #monitor# - 有关此包的更多信息，请参阅 http://runjags.sourceforge.net/quickjags.html

[编辑:监测 T.obs 并不是真的有必要，但这表明 T.obs 中的缺失值都被估计为 > 3 并且观察到的值如预期的那样是非随机的]