# RIGHT-CENSORING EXAMPLE WITH AML DATA # "Survival in patients with Acute Myelogenous Leukemia. The question at the time was # whether the standard course of chemotherapy should be extended ('maintainance') for # additional cycles." status=1 means NOT censored lapply(c("survival","runjags"), library, character.only=TRUE) print(aml) # time status x # 1 9 1 Maintained # 2 13 1 Maintained # 3 13 0 Maintained # 4 18 1 Maintained # 5 23 1 Maintained # 6 28 0 Maintained # 7 31 1 Maintained # 8 34 1 Maintained # 9 45 0 Maintained #10 48 1 Maintained #11 161 0 Maintained #12 5 1 Nonmaintained #13 5 1 Nonmaintained #14 8 1 Nonmaintained #15 8 1 Nonmaintained #16 12 1 Nonmaintained #17 16 0 Nonmaintained #18 23 1 Nonmaintained #19 27 1 Nonmaintained #20 30 1 Nonmaintained #21 33 1 Nonmaintained #22 43 1 Nonmaintained #23 45 1 Nonmaintained # CREATE TRUE/FALSE VECTOR FROM 0/1 VECTOR censored <- !aml$status # JAGS WEIBULL FIT WITH CENSORING # aml.time IS HOW THE DATA ARE RECORD WITH OBSERVED TIMES INCLUDING aml.time <- aml$time # CENSORING TIMES FOR dinterval(t[i], c[i]) SECOND PARAMETER # UNIQUE VERSION REQUIRED FOR FUNCTION aml.cens <- aml$time # aml.na IS TIME VECTOR WITH THE CENSORED TIMES ASSIGNED NA aml.na <- aml$time aml.na[censored] <- NA aml.model <- "model { for(i in 1:PATIENTS) { is.censored[i] ~ dinterval(aml.na[i], aml.cens[i]) aml.na[i] ~ dweib(alpha, beta[group[i]]) surv[i] <- 1 - pweib(aml.time[i], alpha, beta[group[i]]) } alpha ~ dgamma(1, 0.1) beta[1] ~ dgamma(1, 0.1) beta[2] ~ dgamma(1, 0.1) }" aml.data <- dump.format(list( aml.time = aml.time, aml.na = aml.na, aml.cens = aml.cens, is.censored = 1 - aml$status, group = as.numeric(aml$x), PATIENTS = length(aml.na) )) # aml.init IS A VECTOR WITH NA'S FOR UNCENSORED DATA VALUES # AND CENSOR TIME PLUS 1 FOR THE CENSORED DATA VALUES aml.init <- rep(NA, length(aml.na)) aml.init[censored] <- aml.cens[censored] + 1 aml.inits <- c( dump.format(list(alpha = 1.1, beta = c(1,0.001), aml.na = aml.init)), dump.format(list(alpha = 2.3, beta = c(1,0.02), aml.na = aml.init)), dump.format(list(alpha = 5.3, beta = c(1,1.50), aml.na = aml.init)) ) aml.out = run.jags(aml.model, monitor=c("surv", "alpha", "beta","deviance"), data=aml.data, inits=aml.inits, burnin = 50000, sample = 100000, n.chains=3) #Compiling rjags model and adapting for 1000 iterations... #Calling the simulation using the rjags method... # Burning in the model for 50000 iterations... # |**************************************************| 100% # Running the model for 100000 iterations... # |**************************************************| 100% #Simulation complete #Calculating the Gelman-Rubin statistic for 27 variables.... #Note: Unable to calculate the multivariate psrf due to an error calculating #the Gelman-Rubin statistic #The Gelman-Rubin statistic is below 1.05 for all parameters #Finished running the simulation ( aml.summary <- summary(aml.out) ) #Iterations = 51001:151000 #Thinning interval = 1 #Number of chains = 3 #Sample size per chain = 1e+05 # #1. Empirical mean and standard deviation for each variable, # plus standard error of the mean: # # Mean SD Naive SE Time-series SE #alpha 0.9805 0.1807 3.30e-04 0.001975 #beta[1] 0.0258 0.0217 3.96e-05 0.000195 #beta[2] 0.0587 0.0394 7.19e-05 0.000375 #deviance 155.0150 4.6167 8.43e-03 0.041372 #surv[1] 0.8358 0.0710 1.30e-04 0.000536 #surv[2] 0.7775 0.0842 1.54e-04 0.000587 #surv[3] 0.7775 0.0842 1.54e-04 0.000587 #surv[4] 0.7115 0.0964 1.76e-04 0.000624 #surv[5] 0.6519 0.1054 1.92e-04 0.000624 #surv[6] 0.5979 0.1121 2.05e-04 0.000614 #surv[7] 0.5679 0.1153 2.11e-04 0.000607 #surv[8] 0.5396 0.1180 2.15e-04 0.000600 #surv[9] 0.4487 0.1242 2.27e-04 0.000585 #surv[10] 0.4269 0.1251 2.28e-04 0.000585 #surv[11] 0.0861 0.0829 1.51e-04 0.000652 #surv[12] 0.7813 0.0783 1.43e-04 0.000640 #surv[13] 0.7813 0.0783 1.43e-04 0.000640 #surv[14] 0.6826 0.0915 1.67e-04 0.000642 #surv[15] 0.6826 0.0915 1.67e-04 0.000642 #surv[16] 0.5719 0.1003 1.83e-04 0.000553 #surv[17] 0.4804 0.1041 1.90e-04 0.000438 #surv[18] 0.3562 0.1045 1.91e-04 0.000283 #surv[19] 0.3013 0.1025 1.87e-04 0.000251 #surv[20] 0.2663 0.1003 1.83e-04 0.000247 #surv[21] 0.2358 0.0977 1.78e-04 0.000264 #surv[22] 0.1594 0.0871 1.59e-04 0.000318 #surv[23] 0.1478 0.0848 1.55e-04 0.000326 # #2. Quantiles for each variable: # # 2.5% 25% 50% 75% 97.5% #alpha 6.49e-01 0.8540 0.9724 1.0984 1.3536 #beta[1] 3.45e-03 0.0112 0.0197 0.0334 0.0837 #beta[2] 1.22e-02 0.0310 0.0491 0.0757 0.1607 #deviance 1.48e+02 151.6841 154.4606 157.7398 165.5861 #surv[1] 6.70e-01 0.7944 0.8461 0.8879 0.9436 #surv[2] 5.88e-01 0.7259 0.7872 0.8393 0.9132 #surv[3] 5.88e-01 0.7259 0.7872 0.8393 0.9132 #surv[4] 5.02e-01 0.6498 0.7197 0.7816 0.8751 #surv[5] 4.31e-01 0.5823 0.6583 0.7281 0.8381 #surv[6] 3.69e-01 0.5219 0.6022 0.6781 0.8026 #surv[7] 3.37e-01 0.4888 0.5710 0.6499 0.7823 #surv[8] 3.07e-01 0.4578 0.5415 0.6230 0.7624 #surv[9] 2.17e-01 0.3598 0.4463 0.5343 0.6957 #surv[10] 1.97e-01 0.3368 0.4233 0.5125 0.6789 #surv[11] 2.79e-03 0.0255 0.0605 0.1205 0.3092 #surv[12] 6.05e-01 0.7338 0.7903 0.8386 0.9076 #surv[13] 6.05e-01 0.7338 0.7903 0.8386 0.9076 #surv[14] 4.87e-01 0.6236 0.6892 0.7484 0.8418 #surv[15] 4.87e-01 0.6236 0.6892 0.7484 0.8418 #surv[16] 3.69e-01 0.5042 0.5748 0.6426 0.7586 #surv[17] 2.80e-01 0.4077 0.4799 0.5523 0.6839 #surv[18] 1.70e-01 0.2807 0.3505 0.4256 0.5743 #surv[19] 1.26e-01 0.2262 0.2933 0.3677 0.5216 #surv[20] 1.00e-01 0.1923 0.2567 0.3301 0.4867 #surv[21] 7.91e-02 0.1633 0.2246 0.2966 0.4546 #surv[22] 3.45e-02 0.0940 0.1445 0.2093 0.3663 #surv[23] 2.90e-02 0.0841 0.1323 0.1955 0.3514 \end{VM} \section*{AML Survival Plot} \begin{VM} postscript("Class.Multilevel/Images/aml.survival.ps") par(mfrow=c(1,1),mar=c(5,5,2,2),lwd=2,col.axis="white",col.lab="white", col.sub="white", col="white",bg="slategray", cex.lab=1.3) plot(survfit(Surv(time) ~ status, data = aml), lty = 2:3, col = c("red","darkgreen"), lwd=3, xlab = "Time", ylab = "Survival", main = "Acute Myelogenous Leukemia Survival") # POINTS FOR Maintained points(aml$time[as.numeric(aml$x) == 1], aml.summary$statistics[5:27,1][as.numeric(aml$x) == 1], col="pink", pch=19) # POINTS FOR Nonmaintained points(aml$time[as.numeric(aml$x) == 2], aml.summary$statistics[5:27,1][as.numeric(aml$x) == 2], col="lightgreen", pch=19) \end{VM} \section*{AML Survival Plot} \begin{VM} # CREDIBLE INTERVAL FOR Maintained segments(aml$time[as.numeric(aml$x) == 1], aml.summary[1:23,1][as.numeric(aml$x) == 1], aml$time[as.numeric(aml$x) == 1], aml.summary[1:23,3][as.numeric(aml$x) == 1], col="pink",lwd=0.5) # CREDIBLE INTERVAL FOR Nonmaintained segments(aml$time[as.numeric(aml$x) == 2], aml.summary[1:23,1][as.numeric(aml$x) == 2], aml$time[as.numeric(aml$x) == 2], aml.summary[1:23,3][as.numeric(aml$x) == 2], col="lightgreen",lwd=0.5) dev.off() library(rjags, arm, coda) system("echo 'y <- structure(c(' > Class.Multilevel/dogs.model.list") write.table(dogs,"Class.Multilevel/dogs.model.list",sep=",",row.names=FALSE,quote=FALSE, col.names=FALSE, append=TRUE) system("echo '), .Dim = as.integer(c(30, 25)))' >> Class.Multilevel/dogs.model.list") dogs.logit.rjags <- function() { for (j in 1:30){ n.avoid[j,1] <- 0 n.shock[j,1] <- 0 for (t in 2:25){ n.avoid[j,t] <- n.avoid[j,t-1] + (1 - y[j,t-1]) n.shock[j,t] <- n.shock[j,t-1] + y[j,t-1] } for (t in 1:25){ y[j,t] ~ dbin (p[j,t], 1) logit(p[j,t]) <- b.0 + b.1*n.avoid[j,t] + b.2*n.shock[j,t] } } b.0 ~ dnorm (0,0.01) b.1 ~ dnorm (0,0.01) b.2 ~ dnorm (0,0.01) } write.model(dogs.logit.rjags, "Class.Multilevel/dogs.logit.rjags") y <- dogs dogs.model <- jags.model(file="Class.Multilevel/dogs.logit.rjags", data=y,n.chains=1,n.adapt=5000) update(dogs.model,n.iter=500000) dogs.mcmc <- coda.samples(model=dogs.model,variable.names=c("b.0","b.1","b.2"), n.iter=250000) summary(dogs.mcmc) dogs <- read.table("Class.Multilevel/data/dogs2.dat", row.names=1) n.dogs <- nrow(dogs); n.trials <- ncol(dogs) b.0 <- 1.80; b.1 <- -0.35; b.2 <- -0.21 n.sims <- 1000 y.rep <- array (NA, c(n.sims, n.dogs, n.trials)) for (j in 1:n.dogs){ n.avoid.rep <- rep (0, n.sims) n.shock.rep <- rep (0, n.sims) for (t in 1:n.trials){ p.rep <- invlogit (b.0 + b.1*n.avoid.rep + b.2*n.shock.rep) y.rep[,j,t] <- rbinom (n.sims, 1, p.rep) n.avoid.rep <- n.avoid.rep + 1 - y.rep[,j,t] n.shock.rep <- n.shock.rep + y.rep[,j,t] } } #dog <- dogs[,-1] par(mfrow=c(1,2),mar=c(2,2,2,2),oma=c(2,2,4,2),col.axis="white",col.main="white", col.lab="white",col.sub="white",col="white",bg="slategray") plot(1:25,1-apply(dogs,2,mean),type="n") rand.sample <- sample(x=1:n.sims,size=20) for (i in rand.sample) lines(1:25,1-apply(y.rep[i,,],2,mean),col="mintcream") lines(1:25,1-apply(dogs,2,mean),col="darkorchid4",lwd=3,type="l") mtext(side=3,line=1.25,"T(yrep) versus T(y)") plot(1:25,seq(-0.5,0.5,length=25),type="n") for (i in rand.sample) lines(1:25,(1-apply(y.rep[i,,],2,mean))-(1-apply(dogs,2,mean)), col="mintcream",lwd=2,type="l") abline(h=0,col="darkorchid4",lwd=3) mtext(side=3,line=1.25,"T(yrep)−T(y)") title(outer=TRUE,line=2,col="white","Avoidance By Trial Number") dogsort <- function (y) { n.dogs <- nrow(y) n.trials <- ncol(y) last.shock <- rep (NA, n.dogs) for (j in 1:n.dogs) last.shock[j] <- max ((1:n.trials)[y[j,]==1]) print(last.shock) y[order(last.shock),] } dogsort(dogs) y <- dogs dogs.model <- jags.model(file="Class.Multilevel/dogs.log.rjags", data=y,n.chains=1,n.adapt=50000) update(dogs.model,n.iter=500000) dogs.mcmc <- coda.samples(model=dogs.model,variable.names=c("b.1","b.2"), n.iter=250000) summary(dogs.mcmc) dogs <- read.table("http://jgill.wustl.edu/data/dogs2.dat", row.names=1) n.dogs <- nrow(dogs); n.trials <- ncol(dogs) b.1 <- -0.24; b.2 <- -0.08 n.sims <- 1000 y.rep <- array (NA, c(n.sims, n.dogs, n.trials)) for (j in 1:n.dogs){ n.avoid.rep <- rep (0, n.sims) n.shock.rep <- rep (0, n.sims) for (t in 1:n.trials){ p.rep <- exp(b.1*n.avoid.rep + b.2*n.shock.rep) y.rep[,j,t] <- rbinom (n.sims, 1, p.rep) n.avoid.rep <- n.shock.rep + 1 - y.rep[,j,t] n.shock.rep <- n.shock.rep + y.rep[,j,t] } } par(mfrow=c(1,2),mar=c(2,2,2,2),oma=c(2,2,4,2),col.axis="white",col.main="white", col.lab="white",col.sub="white",col="white",bg="slategray") plot(1:25,1-apply(dogs,2,mean),type="n") rand.sample <- sample(x=1:n.sims,size=20) for (i in rand.sample) lines(1:25,1-apply(y.rep[i,,],2,mean),col="mintcream") lines(1:25,1-apply(dogs,2,mean),col="darkorchid4",lwd=3,type="l") mtext(side=3,line=1.25,"T(yrep) versus T(y)") plot(1:25,seq(-0.5,0.5,length=25),type="n") for (i in rand.sample) lines(1:25,(1-apply(y.rep[i,,],2,mean))-(1-apply(dogs,2,mean)), col="mintcream",lwd=2,type="l") abline(h=0,col="darkorchid4",lwd=3) mtext(side=3,line=1.25,"T(yrep)−T(y)") title(outer=TRUE,line=2,col="white","Avoidance By Trial Number, Log Version")