/* SAS code to enter the disease free survival data
   for bone marrow transplant patients */

data set1;
 infile "c:\stat565\data\BMT2.txt";
 input patient group time status agep aged wait fab mtx;
 z1=agep-28;
 z2=aged-28;
 z3=z1*z2;
 if(group=2)then z4=1; else z4=0;
 if(group=3) then z5=1; else z5=0;
 z6=fab;
 z7=12*(wait/365);
 z8=mtx; 
run;

proc print data=set1; run;

proc phreg data=set1;
 model time*status(0)= z1-z8 /ties=efron;
run;



/* Examination of functional form of 
   covariates using marrtingale residulas  */

proc phreg data=set1;
 model time*status(0) = z1-z6 z8/ ties=efron;
 output out=resid_out resmart=mart_res /order=data;

data data_res;
 merge set1 resid_out;

/* gplot has a built in smoother with 
   smoothness ranging from 0 to 100 */

proc gplot data=data_res;
 plot mart_res*z7 /haxis=axis2 vaxis=axis1;
 symbol i=sm60s v=dot h=1.2 w=3;
 axis1 label = (h=2 r=0 a=90 f=swiss "Residuals")
        value = (h=2.0 f=swiss);
 axis2 label = (h=2 f=swiss f=swiss "Waiting Time for Transplant (months)")
        value = (h=2.0 f=swiss);
        label mart_res='Residual';
 title "Martingale Residuals: Waiting Time";
  run;


/* Test proportional hazards criterion 
   using Schoenfeld residuals */

proc phreg data=set1;
    model time*status(0) = z1-z8 /ties=efron; 
    output out=schoenb ressch= schz1-schz8;
	run;

data schoenb; set schoenb; if(time<900); run;

proc gplot data=schoenb;
   plot (schz1-schz8)*time/
       haxis=axis2 vaxis=axis1;
   symbol value=dot i=sm60s h=1.2 w=3;
   axis1 label = (h=2 r=0 a=90 f=swiss )
        value = (h=2.0 f=swiss);
 axis2 label = (h=2 f=swiss)
        value = (h=2.0 f=swiss);
   title "Schoenfeld Residuals";
   run;


/* Examine dfbeta residuals for influence: 
   accounts for change in score residuals 
   when dropping an observation.  Information
   matrix is not adjusted, influence may be
   underestimated.  However, much quicker to
   compute than jackknife residuals (refit 
   Cox model for each subject) */

proc phreg data=set1;
 model time*status(0)= z1-z8/ties=efron;
 output out=dfout dfbeta = zdf1-zdf8;
run;

data dfout;  set dfout;
 ID = _n_; 
 run;

proc gplot data=dfout;
 plot (zdf1-zdf8)*ID /
     haxis=axis2 vaxis=axis1;
 symbol1 v=dot h=1.0;
 axis1 label = (h=2 r=0 a=90 f=swiss )
        value = (h=2.0 f=swiss);
 axis2 label = (h=2 f=swiss 'Patient')
        value = (h=2.0 f=swiss);
 title "Dfbeta values";
run;