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

/* ---------------------------------------------------------
 program: PBC3.SAS
 input:   pbc.dat  DATA A
 does:    Cox regression
        -- plots delta betas and Shoenfeld residuals.
        -- Checks proportionality assumption for all preds
        -- computes Harrell's Z
------------------------------------------------------------ */;

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

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

* automatically prints the date and time on the footnote of output;
FOOTNOTE  "PBC3.SAS      &SYSTIME    &SYSDATE";
TITLE 'Schoenfeld and Delta Beta from Cox Regression';

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

* Designates file name where data is stored;
LIBNAME DATA "d:RA2004\SASwork\";

*********************************************;
DATA A;
  SET DATA.PBC;
 * transform some of the analysis variabels;
  logalb = log(albumin);
  logproth = log(proth);
  logbili =  log(bili);
run;

* Outtest option provides access to the Cov(Beta);
* This is needed to compute the std. delta beta;
proc phreg data=A covout outest=C;
  model time*delta(0)=logbili logalb age
        logproth edema;
  id n;
  output out=B ressch = schbili schalb schage schproth schedema
               dfbeta = dfbili dfalb dfage dfproth dfedema   ;
run;

**************** Delta Betas *****************************;
***  Compute the standard errors of betas;
PROC IML;
USE C;
  READ ALL VAR{ LOGBILI LOGALB AGE
                LOGPROTH EDEMA } INTO Z;
CLOSE C;
EST = (Z[1,]);
COVB = Z[2:6,];
STDERR = (SQRT(VECDIAG(COVB)))`;
PRINT EST COVB STDERR;
        VARNAME = {SE_BILI SE_ALB SE_AGE SE_PROTH SE_EDEMA};
        MATTRIB STDERR COLNAME=VARNAME;
        CREATE STDERR FROM STDERR [COLNAME=VARNAME] ;
        APPEND FROM STDERR [COLNAME=VARNAME];
RUN;
QUIT;

** add stderr estimates to each record;
** use to standardize the delta betas;
data D;
  if _N_ = 1 then set STDERR;
  set B;
  sdfbili = -dfbili / se_bili;
  sdfalb = -dfalb / se_alb;
  sdfage = -dfage / se_age;
  sdfproth = -dfproth / se_proth;
  sdfedema = -dfedema / se_edema;
label   sdfbili = 'Std. Delta Beta'
        sdfalb = 'Std. delta beta'
        sdfage = 'Std. delta beta'
        sdfproth = 'Std. delta beta'
        sdfedema = 'Std. delta beta';
run;


** 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=(axis,legend,pattern,symbol,title,footnote)
  hsize=7.5in vsize=7.5in cback=white htext=3pct
  device=win ftext=swissb;
** for unix, use:  device=kermit;

axis1 label=(angle=270 rotate=90'Std. Delta Beta');
** Plot standardized delta betas;
proc gplot data=D gout=mygraphs ;
  plot sdfbili * logbili/vaxis=axis1 noframe ;
  plot sdfalb * logalb/vaxis=axis1 noframe    ;
  plot sdfage *  age/vaxis=axis1 noframe ;
  plot sdfproth * logproth/vaxis=axis1 noframe ;
  plot sdfedema * n/vaxis=axis1 noframe ;
  plot sdfedema * edema/vaxis=axis1 noframe ;
  symbol1 value=dot h=1.5 pct  i=none;
run;
quit;

******* Schoenfeld Residuals & Harrel's Z *******************;
*** limit the dataset to failure times  ;
title "Calculate Harrel's Z";
data E;  set B;
  if schbili = . then delete;
run;

PROC IML;
USE E;
  READ ALL VAR{ TIME } INTO TIME;
  READ ALL VAR { SCHBILI SCHALB
                SCHAGE SCHPROTH SCHEDEMA } INTO RESIDS;
CLOSE E;
RANKTIME = RANK(TIME);
** Compute correlations (taken from p.112 IML manual);
X=RANKTIME || RESIDS;
N=NROW(X);
SUM=X[+,];
XPX=T(X)*X-T(SUM)*SUM/N;
S=DIAG(1/SQRT(VECDIAG(XPX)));
ALLCORR=S*XPX*S;
** keep corr b/w rank and each var;
CORR=ALLCORR[2:6,1];
HARRELZ = 0.5*LOG((J(5,1,1) + CORR) / (J(5,1,1) - CORR))*SQRT(N-3);
PVAL    = 2 * PROBNORM(-ABS(HARRELZ));
VARNAME = {'log bilirubin', 'log albumin', 'age ', 'log proth', 'edema'};
PRINT "Harrell's Z test for a linear trend in the residuals",
       VARNAME CORR HARRELZ PVAL;
RUN;
QUIT;

goptions reset=(axis,legend,pattern,symbol,title,footnote)
  hsize=7.5in vsize=7.5in cback=white htext=3pct
  device=win ftext=swissb;
  symbol1 value=dot h=1.5 pct  i=sm55s ;

  axis1 label=none;
** Plot Shoefeld Residuals;
proc gplot data=D gout=mygraphs ;
title h=4pct 'log Bilirubin';
  plot schbili*time / vref = 0 vaxis=axis1 noframe;

proc gplot data=D gout=mygraphs ;
title h=4pct 'log Albumin';
  plot schalb*time / vref = 0 vaxis=axis1 noframe;

proc gplot data=D gout=mygraphs ;
title h=4pct 'Age';
  plot schage*time / vref = 0 vaxis=axis1 noframe;

proc gplot data=D gout=mygraphs ;
title h=4pct 'log Prothrombin';
  plot schproth*time / vref = 0 vaxis=axis1 noframe;

proc gplot data=D gout=mygraphs ;
title h=4pct 'Edema';
  plot schedema*time / vref = 0 vaxis=axis1 noframe;

run;
quit;


************************************************;
** plot log(-logS(t)) to check the PH assumption;
************************************************;
goptions reset=(axis,legend,pattern,symbol,title,footnote)
  hsize=7.5in vsize=7.5in cback=white htext=3pct
  device=win ftext=swissb;
symbol1 c=DEFAULT v=NONE i=steplj l=3;
symbol2 c=RED     v=NONE i=steplj l=1;
symbol3 c=GREEN   v=NONE i=steplj l=25;
symbol4 c=BLUE    v=NONE i=steplj l=41;

axis1 label=none;

****************** Bilirubin;
proc phreg data=A;
  model time*delta(0)= logalb age
        logproth edema;
  strata bili (0,1.45,3.25,6.75,100);
  baseline out=B loglogs=lls;
run;
legend1 label=none value=(t=1 "bili<1.45"
        t=2 "1.45<bili<3.25"    t=3 "3.25<bili<6.75"
        t=4 "bili>6.75     ") ;
proc gplot data=B gout=mygraphs;
  plot lls*time = bili / legend = legend1 vaxis=axis1 noframe;
run;
quit;


****************** Albumin ;
proc phreg data=A;
  model time*delta(0)= logbili  age
        logproth edema;
  strata albumin (0,3.1,3.4,3.67,100);
  baseline out=B loglogs=lls;
run;
legend1 label=none value=(t=4 "Alb>3.67"
        t=3 "3.4<Alb<3.67" t=2 "3.1<Alb<3.4"
        t=1 "Alb<3.1");
proc gplot data=B gout=mygraphs;
  plot lls*time = albumin / legend = legend1 vaxis=axis1 noframe;
run;
quit;


****************** Age ;
proc phreg data=A;
  model time*delta(0)= logbili logalb
        logproth edema;
  strata age (0,46,53,61,100);
  baseline out=B loglogs=lls;
run;
legend1 label=none value=(t=4 "Age>61"
        t=3 "53<Age<61"         t=2 "46<Age<53"
        t=1 "Age<46");
proc gplot data=B gout=mygraphs;
  plot lls*time =age /legend = legend1 vaxis=axis1 noframe;
run;
quit;


****************** proth ;
proc phreg data=A;
  model time*delta(0)= logbili logalb age
        edema;
  strata proth (0,10.5,11,11.8,100);
  baseline out=B loglogs=lls;
run;
legend1 label=none value=(t=4 "Proth>11.8"
        t=3 "11<Proth<11.8" t=2 "10.5<Proth<11"
        t=1 "Proth<10.5");
proc gplot data=B gout=mygraphs;
  plot lls*time = proth/legend=legend1 vaxis=axis1 noframe;
run;
quit;


****************** edema;
proc phreg data=A;
  model time*delta(0)= logbili logalb age
        logproth;
  strata edema;
  baseline out=B loglogs=lls;
run;
legend1 label=('Edema') value = (t=3 "1"
        t=2 "0.5" t=1 "0");
proc gplot data=B gout=mygraphs;
  plot lls*time = edema/ legend=legend1 vaxis=axis1 noframe;
run;
quit;




******** Graphs on a single page ******************;
* create titles for the grouped graphs;
GOPTIONS  ROTATE=LANDSCAPE HTITLE=100 PCT
        HSIZE=10in VSIZE=.5 in ftext=swiss nodisplay;
PROC GSLIDE GOUT=mygraphs;
        TITLE h=100pct'Standardized Delta-Beta from Mayo Model';
     RUN;
        TITLE h=100pct'Schoenfeld Resid Plots to Check Proportional Hazards Assump';
     RUN;
        TITLE h=100pct'Log(-log(S)) Plots for Checking Proportionality';
     RUN;


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

filename  Bios180 'd:\RA2004\SASwork\page274.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
*GOPTIONS DEVICE=LXLP4037 HSIZE=11 VSIZE=8.5 ROTATE=LANDSCAPE;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = DATA.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT   2:GPLOT1  3:GPLOT2  4:GPLOT3  5:GPLOT4  6:GPLOT5 7:GSLIDE;
  RUN;
QUIT;


filename  Bios180 'd:/bios180/page280.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
*GOPTIONS DEVICE=LXLP4037 HSIZE=11 VSIZE=8.5 ROTATE=LANDSCAPE;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = DATA.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT6  2:GPLOT7  3:GPLOT8  4:GPLOT9  5:GPLOT10 7:GSLIDE1;
  RUN;
QUIT;


filename  Bios1804 'd:\RA2004\SASwork\page284.eps';
goptions device=psepsf gsfname=Bios1804 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
*GOPTIONS DEVICE=LXLP4037 HSIZE=11 VSIZE=8.5 ROTATE=LANDSCAPE;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = DATA.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT11 2:GPLOT12 3:GPLOT13 4:GPLOT14 5:GPLOT15 7:GSLIDE2;
  RUN;
QUIT;

/*PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = DATA.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT   2:GPLOT1  3:GPLOT2  4:GPLOT3  5:GPLOT4  6:GPLOT5 7:GSLIDE;
  TREPLAY 1:GPLOT6  2:GPLOT7  3:GPLOT8  4:GPLOT9  5:GPLOT10 7:GSLIDE1;
  TREPLAY 1:GPLOT11 2:GPLOT12 3:GPLOT13 4:GPLOT14 5:GPLOT15 7:GSLIDE2;
RUN;
QUIT;*/


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