source('utils.R')

set.seed(12345)

C = 365 # censor length (one year)

conventional = c( 159, 189, 191, 198, 200, 207, 220, 235, 245, 250,
                  256, 261, 256, 261, 265, 266, 280, 343, 356, rep( C, 5 ) )
germ_free = c( 158, 192, 193, 194, 195, 202, 212, 215, 229, 230, 237,
               240, 244, 247, 259, 300, 301, 321, 337, rep( C, 10 ) )

DF_conventional = data.frame( time=conventional )
DF_conventional$deceased = 1
DF_conventional$deceased[20:24] = 0 # patient is still alive

DF_germ_free = data.frame( time=germ_free )
DF_germ_free$deceased = 1
DF_germ_free$deceased[15:24] = 0 # patient is still alive

conventional_sf = make_kaplan_meir_survival_curve(DF_conventional)
germ_free_sf = make_kaplan_meir_survival_curve(DF_germ_free)
#postscript("../../WriteUp/Graphics/Chapter11/ex_11_4_survival_curves.eps", onefile=FALSE, horizontal=FALSE)
plot(conventional_sf, do.points=FALSE,
     xlab='days to death or end of trial',
     ylab='survival probality',
     main='Exercise 11.4',
     col='blue')
plot(germ_free_sf, do.points=FALSE,
     xlab='days to death or end of trial',
     ylab='survival probality',
     main='Exercise 11.4',
     col='purple',
     add=TRUE)
legend( 'topright', c('conventional','germ free'), lty=c(1,1), col=c('blue', 'purple') )
grid()
#dev.off()

# Estimate the confidence interval of the median survival each condition:
#
library(boot)

B = 500
boot_conventional = boot(DF_conventional$time, samplemedian, B)
boot_germ_free = boot(DF_germ_free$time, samplemedian, B)

s = sd(boot_conventional$t)
conv = c( boot_conventional$t0, boot_conventional$t0 - 2*s, boot_conventional$t0 + 2*s )

s = sd(boot_germ_free$t)
gf = c( boot_germ_free$t0, boot_germ_free$t0 - 2*s, boot_germ_free$t0 + 2*s )

print(conv)
print(gf)



# Part (e):
#
print(sprintf('Conventional mean=%5.6f; Germ free mean=%5.6f', mean(DF_conventional$time), mean(DF_germ_free$time)))