last = c(  4,  6, 10, 11, 11, 11, 13, 17,
          20, 20, 21, 22, 24, 24, 29, 30,
          30, 31, 33, 34, 35, 39, 40, 41,
          43, 45, 46, 50, 56, 61, 61, 63,
          68, 82, 85, 88, 89, 90, 93,104,
         110,134,137,160,169,171,173,175,
         184,201,222,235,247,260,284,290,
         291,302,304,341,345);


dead = c(  1,  1,  1,  1,  1,  1,  1,  1,
           1,  1,  1,  1,  1,  1,  1,  1,
           1,  1,  1,  1,  1,  1,  1,  0,
           0,  1,  1,  1,  1,  0,  0,  1,
           1,  1,  1,  1,  1,  1,  1,  1,
           1,  1,  1,  0,  1,  1,  1,  1,
           1,  1,  1,  0,  0,  0,  0,  0,
           0,  0,  0,  0,  0);
       
# homegrown code for Kaplan-Meier (K-M) and 
#                and Nelson-Aalen (N-A) estimators

n.patients=length(last)

event.times = last[dead==1]; n.deaths=length(event.times)
unique.event.times = unique(event.times);  
n.unique.event.times =length(unique.event.times)

c(n.patients, n.deaths, n.unique.event.times)

times.in.lifetable = c(0,unique.event.times);
n.times.in.lifetable = length(times.in.lifetable);

n.at.risk = c(n.patients);
n.who.died =  c(); n.at.risk=c();
for (i in 1 : length(times.in.lifetable) ) {
 n.who.died = c(n.who.died,sum(event.times==times.in.lifetable[i]));
 n.at.risk  = c(n.at.risk, sum(last>=times.in.lifetable[i]));
	}
conditional.surv.probilities = (n.at.risk - n.who.died) / n.at.risk
unconditional.surv.probilities = cumprod(conditional.surv.probilities);
standard.errors = unconditional.surv.probilities * sqrt(
                  cumsum( n.who.died/( n.at.risk*(n.at.risk - n.who.died) ) )
                  );

nelson.aalen = exp( - cumsum(n.who.died/n.at.risk) )
plot(times.in.lifetable,unconditional.surv.probilities,
     type="s",ylab="Surviving Fraction",xlab="Day Post-Diagnosis",
     ylim=c(0,1.05),xlim=c(0,max(last)));
lines(times.in.lifetable,nelson.aalen,type="s",col="red");
lines(times.in.lifetable,
     unconditional.surv.probilities-1.96*standard.errors,type="s",col="gray")
lines(times.in.lifetable,
     unconditional.surv.probilities+1.96*standard.errors,type="s",col="gray")
nelson.aalen = exp( - cumsum(n.who.died/n.at.risk) )
segments( 70,.85,85,.85);            text( 90,.85,"K-M",adj=c(0,0.5))
segments(115,.85,130,.85,col="red"); text(135,.85,"N-A",adj=c(0,0.5),col="red")
                                     text(160,.85,"point estimate",adj=c(0,0.5))

segments(70,.80, 85,.80,col="gray");
  text(90,.80,"K-M point estimate - 1.96SE",adj=c(0,0.5))
segments(70,.90, 85,.90,col="gray");
  text(90,.90,"K-M point estimate + 1.96SE",adj=c(0,0.5))
  
  
# Using 'survival' package

require(survival)

help(survfit)

# K-M

fit=survfit(Surv(last,dead)) ; fit ;
fit$time[1:22] ; fit$surv[1:22] ; 61*fit$surv[1:22]
plot(fit,ylab="Surviving Fraction [K-M]",xlab="Day Post-Diagnosis")

# N-A

fit=survfit(Surv(last,dead),type="fleming-harrington") ; fit ;fit$time[1:22] ; fit$surv[1:22] ; 61*fit$surv[1:22]
plot(fit,ylab="Surviving Fraction[N-A]",xlab="Day Post-Diagnosis")