/* ---------------------------------------------------------
 program: PBC3-v92.SAS
 input:   pbc.dat  DATA A
 does:    Cox regression
        -- plots delta betas and Shoenfeld residuals.
        -- Checks proportionality assumption 
        -- computes Harrell's Z
------------------------------------------------------------ */;

OPTIONS PS=68 LS=120 PAGENO=1 NODATE NOCENTER ERRORS=3;

TITLE 'Schoenfeld and Delta Beta from Cox Regression';

libname out 'H:\bios180\coursepack\SAS-unc-workshop-2008\';
*********************************************;
DATA A;
  SET out.PBC;
 * transform some of the analysis variabels;
  logalb = log(albumin);
  logproth = log(proth);
  logbili =  log(bili);
  trt=(trt1=1);
  label bili="Serum Bilirubin, mg/dl"
        n="case number"
        time="time in days";
run;

* Outtest option provides access to the Cov(Beta);
* This is needed to compute the std. delta beta;
proc phreg data=A covout outest=C;
  model time*delta(0)=logbili logalb age
        logproth edema;
  id n;
  output out=B ressch = schbili schalb schage schproth schedema
               dfbeta = dfbili dfalb dfage dfproth dfedema   ;
run;

**************** Delta Betas *****************************;
***  Compute the standard errors of betas;


PROC IML;
USE C;
  READ ALL VAR{ LOGBILI LOGALB AGE
                LOGPROTH EDEMA } INTO Z;
CLOSE C;
EST = (Z[1,]);
COVB = Z[2:6,];
STDERR = (SQRT(VECDIAG(COVB)))`;
PRINT EST COVB STDERR;
        VARNAME = {SE_BILI SE_ALB SE_AGE SE_PROTH SE_EDEMA};
        MATTRIB STDERR COLNAME=VARNAME;
        CREATE STDERR FROM STDERR [COLNAME=VARNAME];
        APPEND FROM STDERR;
RUN;
QUIT; 


** add stderr estimates to each record;
** use to standardize the delta betas;
data D;
  if _N_ = 1 then set STDERR;
  set B;
  sdfbili = -dfbili / se_bili;
  sdfalb = -dfalb / se_alb;
  sdfage = -dfage / se_age;
  sdfproth = -dfproth / se_proth;
  sdfedema = -dfedema / se_edema;
label   sdfbili = 'Std. Delta Beta'
        sdfalb = 'Std. Delta beta'
        sdfage = 'Std. Delta beta'
        sdfproth = 'Std. Delta beta'
        sdfedema = 'Std. Delta beta';
run;


title h=1.2  'Standardized Delta-Beta from Mayo Model ';
proc sgscatter data=D ;
 plot sdfbili*logbili  sdfalb*logalb sdfage*age	sdfproth*logproth sdfedema*n sdfedema*edema
   /markerattrs=(symbol=circlefilled) ;
run;
quit;

************************************************;
** plot log(-logS(t)) to check the PH assumption;
************************************************;

****************** Bilirubin;
proc phreg data=A;
  model time*delta(0)= logalb age
        logproth edema;
  strata bili (0,1.45,3.25,6.75,100);
  baseline out=B1 loglogs=lls1;
run;



****************** Albumin ;
proc phreg data=A;
  model time*delta(0)= logbili  age
        logproth edema;
  strata albumin (0,3.1,3.4,3.67,100);
  baseline out=B2 loglogs=lls2;
run;

****************** Age ;
proc phreg data=A;
  model time*delta(0)= logbili logalb
        logproth edema;
  strata age (0,46,53,61,100);
  baseline out=B3 loglogs=lls3;
run;


****************** proth ;
proc phreg data=A;
  model time*delta(0)= logbili logalb age
        edema;
  strata proth (0,10.5,11,11.8,100);
  baseline out=B4 loglogs=lls4;
run;


****************** edema;
proc phreg data=A;
  model time*delta(0)= logbili logalb age
        logproth;
  strata edema;
  baseline out=B5 loglogs=lls5;
run;

proc sgplot data=B1;
step y=lls1  x=time /group=bili ;
run;
quit;

proc sgplot data=B2;
step y=lls2  x=time /group=albumin ;
run;
quit;


proc sgplot data=B3;
step y=lls3  x=time /group=age ;
run;
quit;


proc sgplot data=B4;
step y=lls4  x=time /group=proth ;
run;
quit;


proc sgplot data=B5;
step y=lls5  x=time /group=edema ;
run;
quit;


******* Schoenfeld Residuals & Harrel's Z *******************;
*** limit the dataset to failure times  ; 
title "Calculate Harrel's Z";
data E;  set B;
  if schbili = . then delete;
run; 

PROC IML;
USE E;
  READ ALL VAR{ TIME } INTO TIME;
  READ ALL VAR { SCHBILI SCHALB
                SCHAGE SCHPROTH SCHEDEMA } INTO RESIDS;
CLOSE E;
RANKTIME = RANK(TIME);
** Compute correlations (taken from p.112 IML manual);
X=RANKTIME || RESIDS;
N=NROW(X); 
SUM=X[+,];
XPX=T(X)*X-T(SUM)*SUM/N;
S=DIAG(1/SQRT(VECDIAG(XPX)));
ALLCORR=S*XPX*S;
** keep corr b/w rank and each var;
CORR=ALLCORR[2:6,1];
HARRELZ = 0.5*LOG((J(5,1,1) + CORR) / (J(5,1,1) - CORR))*SQRT(N-3);
PVAL    = 2 * PROBNORM(-ABS(HARRELZ));
VARNAME = {'log bilirubin', 'log albumin', 'age ', 'log proth', 'edema'};
PRINT "Harrell's Z test for a linear trend in the residuals",
       VARNAME CORR HARRELZ PVAL;
RUN;
QUIT;

title h=1.2  'Schoenfeld Resid Plots to Check Proportional
  Hazards Assump ';
proc sgscatter data=D ;
 plot (schbili schalb schage schproth schedema)*time
  /markerattrs=(symbol=circlefilled) loess;
run; quit;

*************************************************;
* Use Cumulative sums of Martingale Residuals
* and their transformation to check
* proportional hazards assumption (assess option);
*************************************************;

ods graphics on;
proc phreg data=A;
  model time*delta(0)=logbili logalb age
        logproth edema;
  assess ph /resample seed=19;
run;
ods graphics off;