Step by step right-censored survival analysis in JAGS

This is a sort of follow-up to an earlier post on SE: https://stats.stackexchange.com/questions/70858/right-censored-survival-fit-with-jags

But here, I would like to see a FULL R script (from start to finish) running a survival analysis on right-censored data in JAGS. All the sites I've found require a very high level of proficiency with JAGS so it's hard for me to understand how to get from one line of code to another. I know this is a lot to ask...

Anyway, here are some example survival data. Groups are t1, t2, t3. NAs refer to right-censored data (censor cutoff = 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())

I know this is a lot to ask, but I have spent days trying to piece something together and I keep getting lost/confused. I know that there are now packages to run this sort of analysis, but I really want to learn how to build this from the ground up and on my own! Thank you, readers!

Solution

I don't do a lot of survival analysis (and you don't state which distribution you would like to use for this part - there are different options), but this code should get you started for the interval censoring part:

    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

This uses the runjags package which does some automated convergence etc diagnostics and allows the shorthand use of #data# and #monitor# within the model code rather than the R code - for more info on this package see http://runjags.sourceforge.net/quickjags.html

[edit: It is not really necessary to monitor T.obs but this demonstrates that the missing values in T.obs are all estimated as > 3 and the observed values are non-stochastic as expected]