### Write model files for the LCSS ZIB + PFS Weibull models
### (except Quality of Life, which lacks zero-inflation)
### Format data for WinBUGS and run the models using R2WinBUGS package
### Run univariate models using R2WinBUGS

library(R2WinBUGS)

## Read in data (long format)
long <- read.table("sim_ZIB_data.txt",header=T)

## Create data list for use in WinBUGS
data <- list()
# Number of subjects
data$N <- length(table(long$ID))
# Balanced design, so number of obs per patient
data$J <- max(table(long$ID))
# Total length of data vector
data$NJ <- dim(long)[1]
# Outcome and longitudinal predictors
data$Y <- long$Symptom
data$day <- long$Day
data$treat <- long$Treatment
data$ID <- long$ID
# Outcome and survival predictors
repeated <- seq(from=1,by=data$J,length=data$N)
data$survtrt <- long$Treatment[repeated]
data$SURV <- long$TTE[repeated]
# Censored data have lower bound=censoring time; uncensored have 0
data$SURV.lb <- ifelse(long$Censored[repeated]==1,
	long$Censored[repeated],0)
rm(repeated)

# Prior parameters
data$betamu0 <- data$betamu1 <- rep(0,3)
data$betamu2 <- rep(0,2)
data$Sigma0 <- data$Sigma1 <- diag(3)*.01
data$Sigma2 <- diag(2)*.01
data$R <- diag(2)*2
data$meen <- c(0,0)

# Function to generate initial values
initials <- function(){
   list(SURV = ifelse(is.na(data$SURV),data$SURV.lb+1,NA),
     beta0 = rnorm(3,data$betamu0), 
     beta1 = rnorm(3,data$betamu1), 
     beta2 = rnorm(2,data$betamu2),
     U = matrix(rnorm(data$N*2),data$N,2), 
     TauU = diag(2)*runif(1), alpha=rnorm(2),
     phi = runif(1),gamma = runif(1))
 }

## Fit the model in WinBUGS
fit <- bugs(model.file="model.txt",inits=initials,data=data,
  parameters.to.save=c('beta0','beta1','beta2','phi','gamma','alpha','U'),
  n.chains=3,n.iter=500,n.burnin=100,debug=T)