# homegrown code for Kaplan-Meier (K-M) 

# and illustration of "Redistribute to the Right algorithm (Efron)

# assuming data are sorted smallest->largest, no ties
	
time = c(2, 4, 6, 8 , 10, 12, 14, 16) ; event = c(1,0,1,0,1, 0, 0, 1)

cbind(time, event)

censored.obsns=time[event==0]; censored.obsns
n.censored = length(censored.obsns); n.censored

n.patients=length(time) ; n.patients

event.times = time[event==1]; 

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

n.at.risk =  c();
n.events  =  c(); 

for (i in 1 : n.times.in.lifetable ) {
 n.events   = c(n.events,  sum(event.times==t[i]));
 n.at.risk  = c(n.at.risk, sum(time >= t[i]));
	}

condnl.surv.prob = (n.at.risk - n.events) / n.at.risk
uncondnl.surv.prob = cumprod(condnl.surv.prob);

jump = c(0,  uncondnl.surv.prob [1: (n.times.in.lifetable-1) ] - 
             uncondnl.surv.prob[2:n.times.in.lifetable]  )
cbind(t, n.at.risk, n.events, condnl.surv.prob, uncondnl.surv.prob, jump)


d.y=0.25

plot(t,uncondnl.surv.prob, type="s",
   xlim=c(-0.2, 1.2)*t[n.times.in.lifetable], 
   ylim=c(-d.y*(n.censored+3),1.1) , main="Redistribution-to-Right Algorithm (Efron)")

for (i in 1:n.times.in.lifetable) {
	S = uncondnl.surv.prob[i]
	text(t[i], 
	    S,paste("S=",toString(round(S,3))), adj=c(1.1,0.5),cex=1.25 )
	text(t[i+1], 
	    (uncondnl.surv.prob[i+1] + uncondnl.surv.prob[i])/2,
	  paste("jump=",toString(round(jump[i+1],3))), adj=c(-0.1,0.5) )
	}
	
for (i in 1:n.patients) {
	if(event[i]==1) {points(c(time[i]), c(-0.05), pch=19, cex=1 );
		             points(c(time[i]), c( 1.05), pch=19, cex=1 ) }
		
	if(event[i]!=1) { points(c(time[i]), c(-0.05), pch=21, cex=1 );
		              points(c(time[i]), c( 1.05), pch=21, cex=1 )
		              text(time[i], 1.1, "-->", adj=c(0.5,0),cex=0.75 ) }
	}

		              

#Greenwood SE

standard.errors = uncondnl.surv.prob * sqrt(
                  cumsum( n.events/( n.at.risk*(n.at.risk - n.events) ) )
                  );
standard.errors


# Illustration of "distribute to right" algorithm of Efron
		
	# assuming, as above, that data are sorted smallest->largest, no ties
	
	n=n.patients
	
	# give each a probability mass of 1/n
	
	# successively redistribute mass associated with each censored obsn.
	
	C=0; colour=c("grey50","black")
	
	for ( i in 1:(n-1) ) {
		
		if(i==1) {
			      mass=rep(1/n,n) ;
			      print( "Initially", quote=F )
			      print( cbind(time, event, mass) ) 
			      for (j in 1:n ) text(time[j],-1.5*d.y,toString(round(1/n,3) ),
			                      adj=c(0.5,0), cex=1,col=colour[event[j]+1 ] )
	             } 
		
		if(event[i]==0) {
			 C=C+1
			 print( " ", quote=F )
			 print( paste("Redistribute mass of ", 
			               toString(round(mass[i],3)), 
			               "at time=",
			               toString(time[i])
			              ), 
			        quote=F )
			 print( " ", quote=F )
			 n.sharing = n - i
			 extra= (1/n.sharing)*mass[i]
		     mass[(i+1):n] = mass[(i+1):n] + extra
		      for (j in (i+1):n ) {
			 	text(time[j],-(C+1)*d.y, toString(round(extra,3)),
			 	     cex=0.65,col=colour[event[j]+1 ], adj=c(0.5,0) ) ;
			    
			    text(time[j],-(C+1.5)*d.y, toString(round(mass[j],3)),
			        cex=1,col=colour[event[j]+1 ], adj=c(0.5,0) ) ;
			    }
			  segments(time[i], -(C+0.8)*d.y, time[i], -(C+1)*d.y )
			  segments(time[i], -(C+1)*d.y, 0.9*time[i+1], -(C+1)*d.y )  
		     mass[i]=0
		     result=cbind(time, event, mass)
		     print( result[result[,"mass"] > 0,] )
		     		
			 } # end censored obsn.
			 
    } # i
text( -0.2*t[n.times.in.lifetable], -d.y*(n.censored+3), "jh 2010.09.13" , cex=0.5, adj=c(0,0) )