5

这是对 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())

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

4

1 回答 1

4

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

    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。网络/quickjags.html

[编辑：实际上没有必要监控 T.obs，但这表明 T.obs 中的缺失值全部估计为 > 3，并且观察到的值如预期的那样是非随机的]

于 2015-05-04T06:51:34.830 回答