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

/* ---------------------------------------------------------
 program: PBC2.SAS
 input:   pbc.dat  DATA A
 does:    Plots residuals from Cox regression
     - linear pred * (mart dev resids full model)
     - qq plot of mart and dev resids
     - mart resids full model * each pred
     - null model mart resids * transformed pred
     - null model mart resids * non-transformed pred
     - (log(bili) bili)*(mart dev resids from model w/o bili
     - each pred * mart resids from model w/o that pred
------------------------------------------------------------ */;

*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  "PBC2.SAS      &SYSTIME    &SYSDATE";
TITLE 'Residual Plots';

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

* Designates file name where data is stored;
*LIBNAME DATA "&PATH\DISTRIB\";

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

**************************************************;
*** plot martingale and deviance resids;
*** by the linear predictor (xbeta);
proc phreg data=A ;
  model time*delta(0) = logbili logalb
        age logproth edema;
  output out=B resmart=mart resdev=dev xbeta=linpred;
run;

data B; set B;
  label mart='Martingale Residuals'
        dev ='Deviance Residuals'
        linpred='Linear Predictor';
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)
  cback=white ftext=swissB htext=3pct htitle=6pct
  hsize=4in vsize=8in device=win;
   symbol1 value=dot h=1.5 pct i=none;


axis1 label=(angle=270 rotate=90'Martingale Residuals');
filename  Bios180 'd:/bios180/page253.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
title f=swiss 'Martingale Residuals for Mayo Model';
proc gplot data=B gout=mygraphs ;
  plot mart*linpred / vaxis=axis1 vref=0 noframe;
run;

axis2 label=(angle=270 rotate=90'Deviance Residuals');
filename  Bios180 'd:/bios180/page254.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
title f=swiss 'Deviance Residuals for Mayo Model';
proc gplot data=B gout=mygraphs ;
  plot dev*linpred  / vaxis=axis2 vref=0 noframe;
run;


****************************************************;
*** QQplots of mart and Dev residuals vs. std normal;
axis3 label=(angle=270 rotate=90'Martingale Residuals');
filename  Bios180 'd:/bios180/page255.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
title f=swiss 'Martingale Residuals for Mayo Model';
proc capability data=B graphics noprint gout=mygraphs;
  qqplot mart / vaxis=axis3 normal(mu=0 sigma=1) nolinelegend;
run;

axis4 label=(angle=270 rotate=90'Deviance Residuals');
filename  Bios180 'd:/bios180/page256.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=12 in vsize=10 in ;
title f=swiss 'Deviance Residuals for Mayo Model';
proc capability data=B graphics noprint gout=mygraphs;
  qqplot dev / vaxis=axis4 normal(mu=0 sigma=1) nolinelegend;
run;


**************************************************;
*** get dev and mart resid from a;
*** model that excludes bilirubin;
goptions reset=(axis,legend,pattern,symbol,title,footnote)
  cback=white htext=3pct ftext=swissb
  device=win  hsize=5.25in vsize=4in;

  axis1 label=(angle=270 rotate=90'Martingale Residuals');
proc phreg data=A  noprint;
  model time*delta(0) = logalb
        age logproth edema;
  id bili logbili;
  output out=B resmart=mart resdev=dev;
run;

proc gplot data=B gout=mygraphs;
  plot (mart dev)*(bili logbili)/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;


**************************************************;
*** martingale by each variable in model;
goptions reset=(axis,legend,pattern,symbol,title,footnote)
  hsize=7.5in vsize=7.5in cback=white htext=3pct
  device=win ftext=swissb;

axis1 label=(angle=270 rotate=90'Martingale Residuals');
proc gplot data=B gout=mygraphs ;
  plot mart*(logbili logalb age logproth edema)/vaxis=axis1 noframe ;
  symbol1 value=dot  h=1.5 pct i=sm60s;
run;

axis2 offset=(1cm,1cm);
proc gplot data=B gout=mygraphs;
  plot mart*edema / haxis=axis2 vaxis=axis1 noframe;
  symbol1 value=dot  h=1.5 pct i=sm60s;
run;
   quit;

*************************************************;
****  Martingal Residuals from a Null Model;
****  Trick into fitting null model;
****  by setting max iterations= 0;
proc phreg data=A;
  model time*delta(0) = logbili logalb
        age logproth edema / maxiter=0;
  id  bili albumin proth n;
  output out=B resmart=mart ;
run;

***********   vs. Non-transformed Covariates;
proc gplot data=B gout=mygraphs;
  plot mart*(bili albumin age proth edema)/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

axis2 offset=(1cm,1cm);
proc gplot data=B gout=mygraphs;
  plot mart*edema / haxis=axis2 vaxis=axis1 noframe ;
  symbol1 value=dot  h=1.5 pct i=box;
run;

***********   vs. Transformed Coveriates;
proc gplot data=B gout=mygraphs;
  plot mart*(logbili logalb age logproth edema)/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

axis2 offset=(1cm,1cm);
proc gplot data=B gout=mygraphs;
  plot mart*edema / haxis=axis2 vaxis=axis1 noframe;
  symbol1 value=dot  h=1.5 pct i=box;
run;

************************************************;
****  Martingal Residuals adjusted variable plots;
*** get mart from model that excludes covars;
*** one at a time from the model. Plot ;
*** resids by the excluded covariate. ;

*** logbili;
proc phreg data=A;
  model time*delta(0) =  logalb
                 age logproth edema;
  id logbili ;
  output out=B resmart=mart ;
run;
proc gplot data=B gout=mygraphs;
  plot mart*logbili /vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

*** logalb;
proc phreg data=A;
  model time*delta(0) = logbili
                 age logproth edema;
  id  logalb;
  output out=B resmart=mart ;
run;
proc gplot data=B gout=mygraphs;
  plot mart* logalb/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

*** age;
proc phreg data=A;
  model time*delta(0) = logbili logalb
                  logproth edema;
  id  age ;
  output out=B resmart=mart ;
run;
proc gplot data=B gout=mygraphs;
  plot mart*age/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

*** logproth;
proc phreg data=A;
  model time*delta(0) = logbili logalb
                 age  edema;
  id  logproth;
  output out=B resmart=mart ;
run;
proc gplot data=B gout=mygraphs;
  plot mart*logproth/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

*** edema;
proc phreg data=A;
  model time*delta(0) = logbili logalb
                 age logproth ;
  id  edema ;
  output out=B resmart=mart ;
run;
proc gplot data=B gout=mygraphs;
  plot mart* edema/vaxis=axis1 noframe;
  symbol value=dot  h=1.5 pct i=sm60s;
run;

***************************************;
* create titles for the grouped graphs;


goptions ROTATE=LANDSCAPE 
        htitle=80pct hsize=10in vsize=.3in ftext=swiss nodisplay;

PROC GSLIDE GOUT=mygraphs;
        TITLE h=100pct'Residuals from Mayo Model w/o log(Bilirubin)';
     RUN;
        TITLE h=100pct'Martingale Resids for Mayo Model by each variable';
     RUN;
        TITLE h=100pct'Martingale Resids from Null Model vs. Non-transformed covars ';
     RUN;
        TITLE h=100pct'Martingale Resids from Null Model vs. Transformed covars ';
     RUN;
        TITLE h=100pct'Adjusted Variable plots using Martingale Resids ';
     RUN;


***************************************;
** Combine graphs on single pages       ;
filename  Bios180 'd:/bios180/page258.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=NEW2V;
  TREPLAY 1:GPLOT  2:GPLOT1;
  TREPLAY 1:CAPABILI  2:CAPABIL1;
RUN;
QUIT;
*/
/*filename  Bios180 'd:/bios180/page258.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=append rotate=landscape 
hsize=10.5 in vsize=8.5 in ;*/

PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW4T;
  TREPLAY 1:GPLOT2  2:GPLOT3 3:GPLOT4  4:GPLOT5 5:GSLIDE;
RUN;
QUIT;

filename  Bios180 'd:/bios180/page259.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT6 2:GPLOT7  3:GPLOT8  4:GPLOT9 5:GPLOT10 6:GPLOT11 7:GSLIDE1;
  RUN;
QUIT;

filename  Bios180 'd:/bios180/page261.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT12 2:GPLOT13 3:GPLOT14 4:GPLOT15 5:GPLOT16 6:GPLOT17 7:GSLIDE2;
  RUN;
QUIT;

filename  Bios180 'd:/bios180/page264.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT18 2:GPLOT19 3:GPLOT20 4:GPLOT21 5:GPLOT22 6:GPLOT23 7:GSLIDE3;
  RUN;
QUIT;

filename  Bios180 'd:/bios180/page267.eps';
goptions device=psepsf gsfname=Bios180 gsfmode=replace rotate=landscape 
hsize=10.5 in vsize=8.5 in  display;
PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT24 2:GPLOT25 3:GPLOT26 4:GPLOT27 5:GPLOT28 6:GPLOT29 7:GSLIDE4;
  RUN;
QUIT;
/*PROC GREPLAY NOFS GOUT=GSEG2;
  IGOUT=MYGRAPHS;
  TC = data.TEMPLATE;
  TEMPLATE=NEW6T;
  TREPLAY 1:GPLOT6 2:GPLOT7  3:GPLOT8  4:GPLOT9 5:GPLOT10 6:GPLOT11 7:GSLIDE1;
  TREPLAY 1:GPLOT12 2:GPLOT13 3:GPLOT14 4:GPLOT15 5:GPLOT16 6:GPLOT17 7:GSLIDE2;
  TREPLAY 1:GPLOT18 2:GPLOT19 3:GPLOT20 4:GPLOT21 5:GPLOT22 6:GPLOT23 7:GSLIDE3;
  TREPLAY 1:GPLOT24 2:GPLOT25 3:GPLOT26 4:GPLOT27 5:GPLOT28 6:GPLOT29 7:GSLIDE4;
RUN;
QUIT;*/



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