/******************************************************************
Author: Daniel A. Powers" <dpowers@mail.la.utexas.edu>
		Department of Sociology
		University of Texas at Austin             (512) 471-1122
		Austin, TX 78712-1008                FAX: (512) 471-1748
Date: 29 Sep 1997
This source code is provided free of any charge for public,
 non-commercial use.
This code has not been subject to formal test criteria and is
 provided without any performance guarantees.

Below is a routine to do basic nonparametric survival analysis (no
covariates) using the product-limit or Kaplan-Meier estimator.
I designed it to produce just the essential elements.  You can modify it
to produce other quantities.  I use gnuplot to graph the t vs. S(t), gauss
graphics would be an alternative, but I've never used it.

******************************************************************/

new;
/* example data  duration = t > 0,

                     1  if t is an event time
               d ={
                     0  if t is a right-censored duration
*/

let dta[10,2]=
3.5    1
5.0    0
6.5    1
12.25  1
4.25   0
6.5    1
8.5    0
10.0   0
10.0   1
12.0   1;

t=dta[.,1];d=dta[.,2];

{d1,d2,d3,d4,d5} = pl(t,d);

proc(5)=pl(t,d);
    local dat,c,ut,eps,nrsk,r,i,n,nd,nc,
          nr,st,sest1,sest2,sest,dn,cs,pqn;
dat=t~d;
dat=sortc(dat,1);
t=dat[.,1];d=dat[.,2];c=1-d;
ut=unique((d.==1).*t,1);
eps=1.e-20;
nrsk=rows(dat);
r=zeros(nrsk,1);
   i=1;
     do until i > nrsk;
        if i .== 1; r[1]=nrsk; endif;
        if i .> 1; r[i] = r[i-1] - d[i-1] - c[i-1]; endif;
     i=i+1;
    endo;
n=rows(ut);
nd=zeros(n,1);
nr=zeros(n,1);
st=zeros(n,1);
sest=zeros(n,1);
dn=r.*(r-d);
dn[nrsk]=dn[nrsk-1];
pqn=d./dn;
cs=cumsumc(pqn);

  i=1;
     do until i > n;
         nd[i]=sumc(d.*(t.==ut[i]));
         nr[i]=maxc(r.*(t.==ut[i]));
         nr[1]=nrsk;
         if i==1; st[i]=1; sest[i] = 0; endif;
         if i > 1; st[i]=st[i-1].*(1 -(nd[i]./(nr[i]+1.e-20)));
         sest[i]=st[i].*sqrt( maxc(cs.*(t.==ut[i])) );
       endif;
     i=i+1;
  endo;

format /rd (10,4);
"     Time     Number     Number at  Survival   St. Error  ";
"              Failed       Risk     Function     Surv.    ";
ut~nd~nr~st~sest;
retp(ut,nd,nr,st,sest);
endp;
