/* ---------------------------------------------------------
 program:       PBC1-v92.SAS
 input:         pbc.sas7bdat  
 output:
 does:          Cox regression with stratification
                Dummy vs. True
------------------------------------------------------------ */;

OPTIONS PS=68 LS=80 PAGENO=1 NOCENTER ERRORS=3;

TITLE 'Cox Regression';

libname out 'H:\bios180\coursepack\SAS-unc-workshop-2008\';

**************************************************;
DATA A;
  SET out.PBC;
  trt=(trt1=1);
  label stage = 'Stage of disease';
run;

*********************************************************;
* Now let's compare dummy variable stratification for
* histologic stage of disease to 'true' stratification.  ;

* You get similar results for the RR associated with ;
*       bilirubin level.      ;

* Also check if the baseline hazards in the four;
*       strata appear to be proportional.                 ;
*********************************************************;

title1 'Cox regression for PBC data (Dummy Stratification)';
*Covariate values: male, age of 50, no edema, treatment group,
*bilirubin=3.25, for each stage;
 data D;
  input sex age edema trt bili  stage;
  cards;
  0 50 0 1 3.25 1
  0 50 0 1 3.25 2
  0 50 0 1 3.25 3
  0 50 0 1 3.25 4;
 run;
proc phreg data=A;
  class stage;
  model time*delta(0)=sex age edema
      trt bili stage;
  baseline out=B covariates=D survival=s
      logsurv=ls loglogs=lls cumhaz=chaz;
run;

title2"Estimated Survival Function";
title3"50-yr old male in treatment group, no edema,and bilirubin=3.25 mg/dl";
proc sgplot data=B;
yaxis max=1 min=0;
step y=s x=time /group=stage;
run;
quit;

title2"Estimated Cumulative Hazard Function";
title3"50-yr old male in treatment group, no edema,and bilirubin=3.25 mg/dl";
proc sgplot data=B;
step y=chaz x=time /group=stage;
run;
quit;

*****************************************************************;
title1 "Cox regression for PBC data (true stratification)";

proc phreg data=A;
  model time*delta(0)=sex age edema trt bili;
  strata stage;
  baseline out=C survival=s
    logsurv=ls loglogs=lls cumhaz=chaz;
run; 
     
 title2'Estimated Survival Function';
 title3'at average value of each covariate';
 proc sgplot data=C;
  yaxis max=1 min=0;
  step y=s x=time /group=stage; 
 run; quit;
     
 title2'Estimated Log(Cumulative Hazard)';
 proc sgplot data=C;
  step y=lls x=time /group=stage; 
 run; quit;
     
     
 title2'Estimated Cumulative Hazard Function';
 proc sgplot data=C;
  step y=chaz x=time /group=stage; 
 run; quit;