# S-Plus code to read in breast cancer data, fit reference analyses # with exponential models and compare the two treatments # Also plots and compares with empirical survival functions x<-matrix(scan("breast"),98,3,byrow=T) i<-x[,3]==1 # compute parameters of gamma posterior for the two exponental rates a1<-sum(x[i,2]); b1<-sum(x[i,1]) # treatment 1 = control a2<-sum(x[!i,2]); b2<-sum(x[!i,1]) # treatment 2 = experimental # and look at the means c(a1/b1, a2/b2) # Control treatment looks better # plot the posterior densities t<-seq(0,0.75,length=200) plot(t,dgamma(t*b1,a1)*b1,type='l') lines(t,dgamma(t*b2,a2)*b2,col=3) # sample 20000 from each k<-20000 t1<-rgamma(k,a1)/b1 t2<-rgamma(k,a2)/b2 # difference and ratio hist(t2-t1,prob=T) sum(t2>t1)/k # Control treatment looks better! # How much better is control treatment in terms of survival? # ... overlay the two exponential survival curves # with posterior means for rates "plugged in t<-seq(0,8,length=100) plot(t, exp(-t*a1/b1) , type='l', ylim=c(0,1) ) lines(t, exp(-t*a2/b2) , col=2) # More incisively, look at posterior for relative risks -- # ratio of survival probabilities at x=1, 4 years etc from start of treatment rr <-exp((t2-t1)) hist(rr,prob=T ) sum(rr)/k mean(rr) rr <-exp((t2-t1)*4) hist(rr,prob=T ) mean(rr) # Does exponentail model fit? # Graph "Empirical Survival Curves" and overlay exponentials with # posterior means for rates "plugged in" # Control group plot(survfit( Surv(x[i,1],x[i,2])~x[i,3])) # magic?! lines(t, exp(-t*a1/b1) , col=3) # Experimental group plot(survfit( Surv(x[!i,1],x[!i,2])~x[!i,3])) # more magic! lines(t, exp(-t*a2/b2) , col=3)