/* SAS code to for Kaplan-Meier estimation
   of surviror functions to times from the
   VA lung cancer trial of 137 male patients
   with inoperable lung cancer */

/* Variables
     Treatment: 1=standard, 2=test (chemotherapy)
     Celltype: 1=squamous, 2=smallcell, 3=adeno, 4=large
     Survival in days
     Status: 1=dead, 0=censored
     Karnofsky score:  Measures patient 
      performance of activities of daily 
      living.
                 
     SCORE   FUNCTION
      100    Normal, no evidence of disease
       90    Able to perform normal activity with 
             only minor symptoms
       80    Able to perform normal activity with 
             effort, some symptoms
       70    Able to care for self but unable to do 
             normal activities
       60    Requires occasional assistance
       50    Requires considerable assistance
       40    Disabled, requires special assistance
       30    Severely disabled
       20    Very sick, requires supportive treatment
       10    Moribund

     Months from Diagnosis
     Age in years
     Prior therapy: 0=no, 10=yes
*/
data va;
    infile 'c:\mydocuments\courses\st565\data\va.dat';
    input rx cellt time status karno months age prior_rx;

    prior_rx = prior_rx/10;

    if (cellt=1) then celltype= 'squamous';
    if (cellt=2) then celltype= 'smallcell';
    if (cellt=3) then celltype= 'adeno';
    if (cellt=4) then celltype= 'large';

/* Create dummy variables for cell types */

cell1=0;
cell2=0;
cell3=0;
if (cellt=2) then cell1=1;
if (cellt=3) then cell2=1;
if (cellt=4) then cell3=1;
dummy=1;
run;


/* select covariate values where the baseline 
   survivor function will be estimated */
 
data cov;
 input rx cell1 cell2 cell3 karno months age prior_rx;
 datalines;
 1 1 0 0 80 6 67 1
 0 1 0 0 80 6 67 1
run;


proc phreg data=va;
 model time*status(0)= rx cell1 cell2 cell3 karno 
     months age prior_rx/ties=efron;
   baseline out=ba xbeta=xbeta stdxbeta=stdxbeta 
            survival=surv upper=ucl lower=lcl
            covariates=cov;
run;

proc print data=ba; run;

proc sort data=ba; by rx time;  run;
 
proc gplot data=ba;
 plot surv*time lcl*time ucl*time/overlay;
 symbol1 v=none interpol=step line=1 w=4;
 symbol2 v=none interpol=step line=2 w=4;
 symbol3 v=none interpol=step line=2 w=4;
 by rx;
 title "Estimated Survival Curves";
run;


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

proc phreg data=va;
 model time*status(0) = dummy;
 output out=resid_out resmart=mart_res /order=data;

data data_res;
 merge va resid_out;

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

proc gplot data=data_res;
 plot mart_res*age /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 )
        value = (h=2.0 f=swiss);
 label mart_res='Residual';
 title "Martingale Residulas";
  run;


/* Test proportional hazards criterion 
   using Schoenfeld residuals */

proc phreg data=va;
    model time*status(0) = rx karno age cell1 
             cell2 cell3 months prior_rx /ties=efron; 
    output out=schoenb ressch= schrx schkarno schage;
	run;

proc gplot data=schoenb;
   plot schrx*time schkarno*time schage*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;

/* Use the SCHOEN marco to produce plots
   of scaled schoenfeld residuals and tests */

%include 'c:\mydocuments\courses\st565\sas\therneau\schoen.sas';

/* Include the macro DASPLINE for fitting splines.  
   This macro is used in the SCHOEN macro */

%include 'c:\mydocuments\courses\st565\sas\therneau\daspline.sas';

%schoen(data=va, time=time, event=status, 
        xvars = rx karno age cell1 cell2 cell3 months prior_rx,
		outsch=outs, outbt=scaled, plot=t, points=yes, df=4,
		alpha=.05);


/* 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=va;
 model time*status(0)= rx karno age cell1 
         cell2 cell3 months age prior_rx/ties=efron;
 output out=dfout dfbeta = rx_inf karno_inf age_inf;
run;

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

proc gplot data=dfout;
 plot rx_inf *ID karno_inf*ID age_inf*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;

/* Robust variance estimate using dfbeta 
   residuals: robust to lack of proportionality, 
   incorrect functional form for covariates, 
   and omitted covariates */

/*
%phlev(data=va, time=time, event=status, 
       xvars =rx cell1 cell2 cell3 karno months age
              prior_rx, id=dummy, collapse=Y) 

/* to detect outliers, deviance (or martingale residuals) */
/*
 output out=resid_out  /order=data;

*/