library("survival")
library("ISwR")
data("melanom")

? melanom

## creates an object of class survival 
? Surv

melanom.surv = Surv(time=melanom$days, 
                   event= melanom$status==1, ## TRUE for death
                   type='right')
str(melanom.surv)
melanom.surv

## Kaplan-Meier estimator of survival function
surv.func.km = survfit(melanom.surv~1)

## Disregarding all individuals who survived (right censored data points)
surv.func.brute = survfit(melanom.surv[melanom$status==1]~1)

##
surv.func.vanilla = matrix(nr=length(unique(melanom$days)),nc=2)
surv.func.vanilla[,1] = sort(unique(melanom$days))
for(i in 1:nrow(surv.func.vanilla))
{
  sub = (melanom$status !=1) | (melanom$days > surv.func.vanilla[i,1])
  surv.func.vanilla[i,2] = sum(melanom$days[sub] >= surv.func.vanilla[i,1]) / length(sub)
}

par(lwd=4)
plot(surv.func.km)
lines(surv.func.brute,col=2)
lines(surv.func.vanilla,col=3)
legend(x='bottomleft',col=1:3,lty=rep(1,3),legend=c('Kaplan-Meier','Brute force','Vanilla'))