FILENAME goptions URL "http://192.38.117.59/~linearpredictors/datafiles/goptions.sas";
%include goptions;

FILENAME url URL "http://192.38.117.59/~linearpredictors/datafiles/readPbc3.sas";
%include url;

/*
Calculating means within quintile groups of bilirubin level.
*/

PROC MEANS DATA = pbc3 MEAN NWAY;
    CLASS biligroup;
    VAR bili;
    OUTPUT OUT = means MEAN=biliquintilemean
RUN;

PROC SORT DATA = pbc3;
    BY biligroup;
RUN;

PROC SORT DATA = means;
    BY biligroup;
RUN;

DATA pbc3;
    MERGE pbc3 means;
    BY biligroup;
    ARRAY q {*} bq1 - bq5; /* dummy variables */
    DO i = 0 TO 4;
        q(i+1) = (biligroup=i);
        END;
RUN;

/*
Created 0-1 variables according to quintile groups of bilirubin.
The next part is calculating Hazard Rate Ratio.
*/

PROC PHREG DATA = pbc3;
    MODEL followup*status(0) = bq2 - bq5;
    OUTPUT OUT = estimates XBETA = loghazardratio / ORDER = DATA;
RUN;

/*
Notice that we requested the estimates data set ordered as pbc3,
so that we do not need a BY statement. 
*/

DATA pbc3;
    MERGE pbc3 estimates;
    KEEP loghazardratio biliquintilemean;
RUN;

/*
Collapsing dataset containing estimates, so that each point only occurs once.
We keep the first observation for each value of biliquintilemean. 
*/

DATA estimates;
  SET pbc3;
  KEEP biliquintilemean loghazardratio;
  BY biliquintilemean;
  IF first.biliquintilemean;
RUN;

/*
The final plot. 
*/

PROC GPLOT UNIFORM DATA = estimates;
    PLOT loghazardratio*biliquintilemean /
        HAXIS = AXIS1 VAXIS = AXIS2; 
    AXIS1
        MINOR = NONE
        ORDER = (0 TO 160 BY 20);      
    AXIS2
        MINOR = NONE
        ORDER = (-1 TO 3 BY 0.5)
        LABEL = (A=90 R=0 'log(Hazard rate ratio)');
    SYMBOL1 VALUE = CIRCLE COLOR = BLACK;
RUN;
QUIT;