################################################
##
## Calculations for Life Table Estimate Example
##
################################################

n.t <- c(80,66,54,39)
d.t <- c(12,8,10,5)
w.t <- c(2,4,5,4)

my.life.table <-data.frame(n.t, d.t, w.t)

# Approach 1:
#m.t <- d.t/n.t

# Approach 2:
#m.t <- d.t/(n.t-w.t)

# Approach 3:
m.t <- d.t/(n.t-w.t/2)

m.t.compl <- 1-m.t

S.hat.t <- cumprod(1-m.t)

my.full.life.table <- cbind(my.life.table, m.t, m.t.compl, S.hat.t)
print(my.full.life.table)

# Using Greenwood's Formula:

sum.elements <- d.t/((n.t-w.t/2)*(n.t-d.t-w.t/2))
var.S.hat.t <- (S.hat.t^2)*cumsum(sum.elements)
se.S.hat.t <- sqrt(var.S.hat.t)

my.full.life.table.se <- cbind(my.full.life.table, se.S.hat.t)
print(my.full.life.table.se)


###################################################
##
## Calculations for Kaplan-Meier Estimate Examples
##
###################################################

library(survival)  # loading the 'survival' package

# Simple Example:

time <- c(2.5,5.5,6.5,9.5,11.5,13.5)

cens.status <- c(1,1,0,1,0,1)  # 1=observed survival time, 0=censoring time

my.surv.fit <- survfit(Surv(time, cens.status) ~ 1, type="kaplan-meier", error="greenwood", conf.type="plain")

summary(my.surv.fit)

plot(my.surv.fit)


### Example 2 (lung cancer patients)

library(survival)

surv.time <- cancer$time
cens.status <- cancer$status - 1   # Converts given (1,2) variable to (0,1)

my.surv.fit.lc <- survfit(Surv(surv.time, cens.status) ~ 1, type="kaplan-meier", error="greenwood", conf.type="plain")

summary(my.surv.fit.lc)

plot(my.surv.fit.lc)


### Example 3 (heart transplant patients)

library(survival)

surv.time <- stanford2$time
cens.status <- stanford2$status   # variable already in (0,1) status

my.surv.fit.ht <- survfit(Surv(surv.time, cens.status) ~ 1, type="kaplan-meier", error="greenwood", conf.type="plain")

summary(my.surv.fit.ht)

plot(my.surv.fit.ht)

# Note:  Other confidence levels could be specified, e.g.:
# my.surv.fit.ht <- survfit(Surv(surv.time, cens.status) ~ 1, type="kaplan-meier", error="greenwood", conf.type="plain", conf.int=0.90)

# Showing point and interval estimates for the median survival time:

my.surv.fit.ht


###################################################
##
## Calculations for log-rank test Examples
##
###################################################

### Example 1 (Ovarian cancer data)

library(survival)

surv.time <- ovarian$futime
cens.status <- ovarian$fustat  ## 1=death and 0=censored
treat.ind <- ovarian$rx   ### 1=cyclophosphamide alone, 2=cyclophosphamide + adriamycin

## Finding the K-M estimator separately for each treatment group:
my.surv.fit.ov <- survfit(Surv(surv.time, cens.status) ~ treat.ind, type="kaplan-meier", error="greenwood", conf.type="plain")
my.surv.fit.ov
plot(my.surv.fit.ov)

## Testing whether treatment 1 is different than treatment 2
## in terms of survival functions

survdiff(Surv(surv.time, cens.status) ~ treat.ind)


### Example 2 (Colon cancer data)

library(survival)

colon.short <- colon[(colon$etype==1) & (colon$rx !='Obs'),]

# colon.short

surv.time <- colon.short$time
cens.status <- colon.short$status  ## 1=recurrence and 0=censored
treat.ind <- ifelse(colon.short$rx=='Lev',1,2)   ### 1="Levamisol", 2="Levamisol + 5FU"

## Finding the K-M estimator separately for each treatment group:
my.surv.fit.col <- survfit(Surv(surv.time, cens.status) ~ treat.ind, type="kaplan-meier", error="greenwood", conf.type="plain")
my.surv.fit.col
plot(my.surv.fit.col)

## Testing whether treatment 2 is better than treatment 1
## in terms of survival functions

survdiff(Surv(surv.time, cens.status) ~ treat.ind)