# 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)