/* PubH 7450
Exampel 8.2 (continued): Estimating the baseline survival and other 
          model-based survival curves from a fitted PHM using the 
          larynx cancer data; see p.284. 
*/

options ls=80 center pagesize=60 number label;
goptions device = PS;  
         /* graphs are saved in a file names sasgraph.ps */

data larynx;
infile '/home/merganser/weip/public_html/course/7450/data/larynx.txt'
       firstobs=19;
input stage time age year status ;
  if stage=1 then Z1=1;
     else Z1=0;
  if stage=2 then Z2=1;
     else Z2=0;
  if stage=3 then Z3=1;
     else Z3=0;

/* covariate values at which survival curves are to be estimated: */
data covvalues;
  input Z1 Z2 Z3 age;
  datalines;
1 0 0 60
0 1 0 60
0 0 1 60
0 0 0 60
;

/* nomean: no estimated survival curve for each covariate at its mean,
           which is given by default. 
   out: save the estimated survival curves in the named SAS dataset
*/
proc phreg data=larynx;
model time*status(0) = Z1 Z2 Z3 age ;
baseline covariates=covvalues out =Pred1 survival=S lower=CIlow upper=CIup
         / nomean;

/* you can then print out the estimated survival curves: */
proc print data=Pred1;
run;

/* draw the estimated survival curves:*/
 
   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(4, .8) 
      value=(f=swiss h=.8 'Stage I' 'Stage II' 
             'Stage III' 'Stage IV'); 
   axis1 label=(h=1 f=swiss a=90) minor=(n=1); 
   axis2 label=(h=1 f=swiss 'Survival Time in Years') minor=(n=4);

 
   proc gplot data=Pred1; 
      plot S*Time=Pattern / legend=legend1 vaxis=axis1  
                            haxis=axis2    cframe=ligr; 
      symbol1 interpol=stepLJ h=1 v=square  c=blue; 
      symbol2 interpol=stepLJ h=1 v=diamond c=yellow; 
      symbol3 interpol=stepLJ h=1 v=circle  c=red; 
      symbol4 interpol=stepLJ h=1 v=dot  c=black; 
      note f=swiss h=1.5 j=c 'Survival of a 60-yr larynx cancer patient'; 
   run;
*********************************************
     Obs    Z1    Z2    Z3    age    time       S        CIlow       CIup

       1     1     0     0     60     0.0    1.00000     .          .     
       2     1     0     0     60     0.1    0.99485    0.98460    1.00000
       3     1     0     0     60     0.2    0.98963    0.97447    1.00000
       4     1     0     0     60     0.3    0.97358    0.94743    1.00000
       5     1     0     0     60     0.4    0.96803    0.93847    0.99851
       6     1     0     0     60     0.5    0.96235    0.92940    0.99647
       7     1     0     0     60     0.6    0.95666    0.92044    0.99431
       8     1     0     0     60     0.7    0.95093    0.91153    0.99204
......
      35     1     0     0     60     7.8    0.46588    0.30084    0.72148
      36     0     1     0     60     0.0    1.00000     .          .     
      37     0     1     0     60     0.1    0.99409    0.98202    1.00000
      38     0     1     0     60     0.2    0.98810    0.96990    1.00000
......