options nocenter nodate ls=80;

** READ DATA **;
data indat;
infile 'c:\zucker\survival\mrfit.dat';
input age 1-2 dbp 3-5 sbp 6-8 race 9 diab 10 chol 11-13 cigs 14-15
  exyr 16-17 exmo 18-19 exday 20-21 fyrs 22-23 icd9 24-27 
  dyr 28-29 dmo 30-31 dday 32-33 died 34;
exdate = mdy(exmo,exday,exyr);
ddate = mdy(dmo,dday,dyr);
fdate = mdy(12,31,90);
time = min(fdate,ddate) - exdate; 
time = time / 365.25;
diab = 2-diab;
white = (race=1);
map = (1/3)*(sbp - dbp) + dbp;
smoke = (cigs gt 0);
if map=. then delete;
died10 = ((died=1) and (time le 10));

** RUN LOGISTIC MODEL AND PUT OUT STATISTICS **;
proc logistic descending;
model died10 = age map diab chol smoke / covb rsq clodds=wald;
units age=5 map=5 chol=10;
output 
  out=odat 
  p=predval 
  c=cval
  cbar=cbval  
  dfbetas=dfbint dfbage dfbmap dfbdiab dfbchol dfbsmok
  difchisq=dchi
  difdev=ddev
  h=hval;
run;

** SIMPLE STATS ON THE DIAGNOSTICS **;
proc means;
  var cval cbval dfbint dfbage dfbmap dfbdiab dfbchol dfbsmok 
    dchi ddev hval;
run;

** NORMALIZE THE DIAGNOSTICS **;
data next;
set; 
cval     =     (cval      -      0.0011657)  /  0.0023549;
cbval    =     (cbval     -      0.0011628)  /  0.0023340;
dfbint   =     (dfbint    -   -1.938665E-7)  /  0.0138895;
dfbage   =     (dfbage    -   -2.000123E-7)  /  0.0138895;
dfbmap   =     (dfbmap    -   -3.543001E-7)  /  0.0141292;
dfbdiab  =     (dfbdiab   -   4.6824824E-8)  /  0.0142767;
dfbchol  =     (dfbchol   -   1.6442058E-6)  /  0.0136560;
dfbsmok  =     (dfbsmok   -   -1.678782E-7)  /  0.0141463;
dchi     =     (dchi      -      1.0024381)  /  1.5726473;
ddev     =     (ddev      -      1.0134396)  /  1.0823136;
hval     =     (hval      -      0.0011190)  /  0.0011400;

** UNIVARIATE STATISTICS ON NORMALIZED DFBETA VALUES **;
proc univariate plot;
  var cval cbval dfbint dfbage dfbmap dfbdiab dfbchol dfbsmok dchi ddev hval;
run;

** RESIDUAL PLOT: NORMALIZED GROUPED RESIDUALS VS. MAP **;
proc rank; var map; ranks maprnk;
data groups; set; 
rnk = 50 * (maprnk / 5362); 
* In this example, 5362 is the sample size, and 50 is the number of groups
* the sample was split up into for the purpose of doing the plots. These
* numbers should be changed to fit the particular application.;
grp = int(rnk-.001) + 1;
proc sort; by grp;
proc means noprint; by grp;
var map predval died10 map;
output out=gmeans mean = zmean pmean ymean n(map)=ng;
data smeans; set gmeans;
expnum = ng * pmean;
obsnum = ng * ymean;
resid = (ymean - pmean);
pvar = pmean * (1-pmean) / ng;
nresid = resid / sqrt(pvar);
lym = log(ymean/(1-ymean));
lpm = log(pmean/(1-pmean));
lresid = lym - lpm;
drop pvar lym lpm;

** GRAPH OF *NORMALIZED* RESIDUALS VERSUS MAP **;
** THIS IS TO CHECK FOR OUTLIER RESIDUALS     **;
symbol1 color=black value=dot;
proc gplot;
plot nresid*zmean;
run;

** GRAPH OF *UNNORMALIZED* LOGISTIC RESIDUALS VERSUS MAP **;
** THIS IS TO CHECK FOR NONLINEAR TREN         D         **;
symbol1 color=black value=dot;
proc gplot;
plot lresid*zmean;
run;

** CREATE NONPARAMETRIC REGRESSION CURVE BASED ON ABOVE PLOT **;
proc loess;
model lresid = zmean / degree=2 direct scale=sd 
  smooth=0.2 0.4 0.6 0.8 1.0;
ods output outputstatistics=ldat1;
run;

** PLOT THE ORIGINAL POINTS AND THE SMOOTH CURVE ON THE SAME GRAPH **;
proc sort data=ldat1;
by smoothingparameter zmean;
symbol1 color=black value=dot;
symbol2 color=black interpol=spline value=none;
proc gplot data=ldat1;
by smoothingparameter;
plot (depvar pred) * zmean / overlay;
run;