** Clears the output and log windows and recalls the program;
DM 'OUT;CLEAR;LOG;CLEAR;PGM;RECALL';

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

*Creates a macro var for location of project work;
%LET PATH=A:        ;

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

* automatically prints the date and time on the footnote of output;
TITLE 'Cox Regression';

* Designates file names to store the program, the output and the log;
FILENAME PROGRAM "d:\bios180\PBC1.SAS";
FILENAME LOG     "d:\bios180\PBC1.LOG";
FILENAME OUTPUT  "d:\bios180\PBC1.OUT";
* Saves the program to above designated file;
DM 'PGM; FILE PROGRAM REPLACE';

* Designates file name where data is stored;
*LIBNAME data "d:\bios180\";

**************************************************;
DATA A;
  SET out.PBC;
  * create dummy vars for stage levels;
  stage2=(stage=2);
  stage3=(stage=3);
  stage4=(stage=4);
  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 stage)';
proc phreg data=A;
  model time*delta(0)=sex age edema
        trt bili stage2 stage3 stage4;
  baseline out=B survival=s
                logsurv=ls loglogs=lls;
run;

****************************************;
*********** Set graphic control information;
* delete work graphics catalog from previous runs;
proc greplay nofs  igout=mygraphs;
   delete _all_; quit;
proc greplay nofs igout=gseg2;
   delete _all_; quit;
  goptions reset=(symbol, pattern, footnote) graphrc
        hpos=0 vpos=0  htext=3 pct ftext=swissb
        ctext= target= gaccess= gsfmode= cback=white;
* goptions device=WIN norotate;  ***for printing to the screen;
  symbol1 c=blue value=none interpol=steplj l=1;
  axis1 label=none;
  axis2   order=0 to 5000 by 1000;
  axis3  label=none order=0 to 1.0  by 0.2;
************;

filename  Bios180 'd:/bios180/page235.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
proc gplot data=B gout=mygraphs;
title1  h=6 pct "Cox regression for PBC data (Dummy stage)";
title2 'Estimated Survival Function';
  plot s*time / haxis=axis2 vaxis=axis3 noframe;
run;


filename  Bios180 'd:/bios180/page236.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;

title2 'Estimated Cumulative Hazard';
data C;
  set B;
 ls=-ls;
run;
proc gplot data=C  gout=mygraphs;
  plot ls*time /  haxis=axis2 vaxis=axis1 noframe;
run;

quit;


*****************************************************************;
title1 "Cox regression for PBC data (true stage)";
proc phreg data=A;
  model time*delta(0)=sex age
        edema trt bili;
  strata stage;
  baseline out=B survival=s
                logsurv=ls loglogs=lls;
run;

*********** Set graphic control information;
goptions reset=( symbol);
symbol1 c=DEFAULT v=CIRCLE   i=steplj l=1;
symbol2 c=RED     v=TRIANGLE i=steplj l=1;
symbol3 c=GREEN   v=SQUARE   i=steplj l=1;
symbol4 c=BLUE    v=STAR     i=steplj l=1;
********** ;

filename  Bios180 'd:/bios180/page237.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
proc gplot data=B gout=mygraphs;
title1  h=6 pct "Cox regression for PBC data (true stage)";
title2 'Estimated Survival Function';
  plot s*time=stage /  haxis=axis2 vaxis=axis3 noframe ;
run;


title2 'Estimated Cumulative Hazard';
data C;
  set B;
 ls=-ls;
run;
proc gplot data=C  gout=mygraphs;
  plot ls*time=stage /  haxis=axis2 vaxis=axis1;
run;

axis1 label=none;
filename  Bios180 'd:/bios180/page238.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
title2 'Log(Cumulative Hazard)';
proc gplot data=B  gout=mygraphs;
  plot lls*time=stage /  haxis=axis2 vaxis=axis1 noframe ;
run;

quit;


***************************************;
** Combine graphs on single pages       ;

GOPTIONS DEVICE=win ROTATE=LANDSCAPE
        HSIZE=10.67in VSIZE=8in;
*GOPTIONS DEVICE=LXLP4037 ROTATE=LANDSCAPE
        HSIZE=11 VSIZE=8.5;*Print in student room;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = DATA.TEMPLATE;
  TEMPLATE=NEW2H;
  TREPLAY 1:GPLOT  2:GPLOT1;
  TREPLAY 1:GPLOT2  2:GPLOT4;
RUN;
QUIT;

** Automatically saves contents of the output and log windows;
DM 'LOG; FILE LOG    REPLACE';
DM 'OUT; FILE OUTPUT REPLACE';