* Read in the data. ;

PROC IMPORT DATAFILE = 'C:\Documents and Settings\richc\My Documents\CPH931F08\LARYNXdata.xls'
     OUT = LD  DBMS = EXCEL REPLACE;
     SHEET =  LARYNXdata;
     GETNAMES = YES;
RUN;

PROC PRINT DATA = LD;
RUN;


* Create indicators for the first three categories of stage. ;

DATA LD;
SET LD;
Z1 = 0;
Z2 = 0;
Z3 = 0;
IF STAGE = 1 THEN Z1 = 1;
IF STAGE = 2 THEN Z2 = 1;
IF STAGE = 3 THEN Z3 = 1;
RUN;


* Begin transcript. ;

ODS RTF FILE = 'C:\Documents and Settings\richc\My Documents\CPH931F08\LARYNXExamples.rtf';

* Fit a Weibull model with log of time to death as the response
* variable and stage as a categorical explanatory variable, in
* the presence of right censoring. ; 

 proc lifereg data=ld;
      class stage;
      model timetodeath*diedinstudy(0)=stage / distribution = weibull;
      output out=new cdf=prob;
   run;

* Plot estimated survival functions for the four strata determined
* by the categorical explanatory variable of stage. ;

   legend1 frame cframe=ligr cborder=black 
      position=center value=(justify=center);

   axis1 label=(angle=90 rotate=0 'Estimated CDF') minor=none;
   axis2 minor=none;

   symbol1 h=1 v=circle c=green i=spline;
   symbol2 h=1 v=circle c=yellow i=spline;
   symbol3 h=1 v=circle c=red i=spline;
   symbol4 h=1 v=circle c=blue i=spline;

   proc sort data=new;
      by prob;

   proc gplot data=new;
      plot prob*timetodeath=stage/ frame cframe=ligr 
          legend=legend1 vaxis=axis1 haxis=axis2;	  
		title 'Weibull Model';
   run;


* Fit a log logistic model with log of time to death as the response
* variable and stage as a categorical explanatory variable, in
* the presence of right censoring. ; 

 proc lifereg data=ld;
      class stage;
      model timetodeath*diedinstudy(0)=stage / distribution = llogistic;
      output out=new cdf=prob;
   run;

* Plot estimated survival functions for the four strata determined
* by the categorical explanatory variable of stage. ;

   legend1 frame cframe=ligr cborder=black 
      position=center value=(justify=center);

   axis1 label=(angle=90 rotate=0 'Estimated CDF') minor=none;
   axis2 minor=none;

   symbol1 h=1 v=circle c=green i=spline;
   symbol2 h=1 v=circle c=yellow i=spline;
   symbol3 h=1 v=circle c=red i=spline;
   symbol4 h=1 v=circle c=blue i=spline;

   proc sort data=new;
      by prob;

   proc gplot data=new;
      plot prob*timetodeath=stage/ frame cframe=ligr 
          legend=legend1 vaxis=axis1 haxis=axis2;
	  title 'Log Logistic Model';
   run;


* Define strata for which estimated survival functions
* are desired with PROC PHREG. ;

   data Inrisks;
      input Z1 Z2 Z3;
      datalines;
   1 0 0
   0 1 0
   0 0 1
   0 0 0
   ;


* Fit an ordinary proportional hazards regression model
* with time to death as the response variable and indicators
* for the first three categories of stage as dichotomous
* explanatory variables, in the presence of right censoring. ;

   proc phreg data=ld  ;
      model timetodeath*diedinstudy(0)=Z1 Z2 Z3;
      baseline covariates=Inrisks out=Pred1 survival=S / nomean;
   run; 


* Plot estimated survival functions for the four strata 
* defined earlier. ;

   data Pred1;
      set Pred1;
      if       Z1 = 1  then Pattern = 1;
      else if  Z2 = 1  then Pattern = 2;
      else if  Z3 = 1  then Pattern = 3;
      else  Pattern = 4;

   legend1 label=none shape=symbol(3, .8)
      value=(f=swiss h=.8 'Stage 1'
             'Stage 2' 'Stage 3' 'Stage 4');
   axis1 label=(h=1 f=swiss a=90) minor=(n=1);
   axis2 label=(h=1 f=swiss) minor=(n=4);
 
   data Pred1;
   set Pred1;
   EstimatedCDF = 1 - S;
   run;

   proc gplot data=Pred1;
      plot EstimatedCDF*TimetoDeath=Pattern / legend=legend1 vaxis=axis1 
                            haxis=axis2    cframe=ligr;
      symbol1 interpol=stepLJ h=1 v=circle  c=green;
      symbol2 interpol=stepLJ h=1 v=circle c=yellow;
      symbol3 interpol=stepLJ h=1 v=circle  c=red;
      symbol4 interpol=stepLJ h=1 v=circle  c=blue;
      title 'Proportional Hazards Model';
   run;

* End transcript. ;

ODS RTF CLOSE;
RUN;