/* PubH 7450
Exampel 11.1b: Diagostics for Cox PHM for the bone marrow transplant data; 
*/

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

/*
g--Disease Group
1-ALL
2-AML Low Risk
3-AML High Risk

T1 -- Time To Death Or On Study Time

T2 --Disease Free Survival Time (Time To Relapse, Death Or End Of Study)

d1 -- Death Indicator
1-Dead 0-Alive

d2 -- Relapse Indicator
1-Relapsed, 0-Disease Free

d3--Disease Free Survival Indicator
1-Dead Or Relapsed, 0-Alive Disease Free)

TA--Time To Acute Graft-Versus-Host Disease

dA--Acute GVHD Indicator
1-Developed Acute GVHD 0-Never Developed Acute GVHD)

TC--Time To Chronic Graft-Versus-Host Disease

dC--Chronic GVHD Indicator
1-Developed Chronic GVHD 0-Never Developed Chronic GVHD

TP--Time To Return of Platelets to Normal Levels

dP--Platelet Recovery Indicator
1-Platelets Returned To Normal, 0-Platelets Never Returned to Normal

Z1--Patient Age In Years

Z2--Donor Age In Years

Z3--Patient Sex
1-Male, 0-Female

Z4--Donor Sex
1-Male, 0-Female

Z5--Patient CMV Status
1-CMV Positive, 0-CMV Negative

Z6--Donor CMV Status
1-CMV Positive, 0-CMV Negative

 Z7--Waiting Time to Transplant In Days

Z8--FAB
1-FAB Grade 4 Or 5 and AML, 0-Otherwise

Z9--Hospital
1-The Ohio State University, 2-Alferd , 3-St. Vincent, 4-Hahnemann

Z10--MTX Used as a Graft-Versus-Host- Prophylactic
1-Yes 0-No
*/

data bmt;
infile '/home/merganser/weip/public_html/course/7450/data/bmt.txt'
       firstobs=87;
input grp T1 T2 d1 d2 d3 Ta da Tc dc Tp dp V1-V10;
if grp=2 then Z1=1; /* AML Low Risk*/
  else Z1=0;
if grp=3 then Z2=1; /* AML High Risk*/
  else Z2=0;
Z3=V8; /* FAB */
Z4=V1-28; /* patient age */
Z5=V2-28; /* donor age */
Z6=Z4*Z5;

/* Output various residuals and some diagnostic stats. 
*/
proc phreg data=bmt;
model T2*d3(0) = Z1 Z2 Z3 Z4 Z5 Z6 ;
output out=diag DFBETA=DBz1 DBz2 DBz3 DBz4 DBz5 DBz6  
       RESSCH=rSCHz1 rSCHz2 rSCHz3 rSCHz4 rSCHz5 rSCHz6 
       RESSCO=rSCOz1 rSCOz2 rSCOz3 rSCOz4 rSCOz5 rSCOz6 
       WTRESSCH=rWTz1 rWTz2 rWTz3 rWTz4 rWTz5 rWTz6 
       LD=ld RESMART=rMART RESDEV=rDEV XBETA=riskscore;

/*
proc print data=diag;
*/
data diag;
  set diag;
  Obs = _N_;

   axis1 label=(h=1 f=swiss a=90) minor=(n=1);

proc gplot data=diag;
  plot rMART*Z5 DBz1*Obs ld*Obs rDEV*riskscore rMART*riskscore / vaxis=axis1;

run;