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;

/* create new variables relating to normal range: 0-17.1 */

DATA pbc3; SET pbc3;
   logbili=log(bili); 
   lbilicent=logbili-log(17.1);
   bililarge=(bili>=17.1); bililargepred=bililarge*tment;
   logbilipos=bililarge*lbilicent; logbilipospred=logbilipos*tment;
RUN;

PROC PHREG DATA=pbc3 OUTEST=esti COVOUT;
    MODEL followup*status(0)=tment bililarge logbilipos
                             bililargepred logbilipospred;        
RUN;

DATA esti; SET esti;
     IF (_NAME_='bililarge' OR _NAME_='logbilipos') THEN DELETE;
     IF _NAME_='followup' THEN DO;
        btment=tment; blp=bililargepred; bpp=logbilipospred; END;
     RETAIN btment blp bpp;        
     IF _NAME_='tment' THEN DO;
        vt=tment; ctl=bililargepred; ctp=logbilipospred; END;
     RETAIN vt ctl ctp;
     IF _NAME_='bililargepred' THEN DO;
        vl=bililargepred; clp=logbilipospred; END;
     RETAIN vl clp;     
     IF _NAME_='logbilipospred' THEN DO;
        vp=logbilipospred; END;
     RETAIN vp;
     OUTPUT;
RUN;

DATA esti; SET esti;
    WHERE _NAME_='logbilipospred';
    key=1;
RUN;

PROC PHREG DATA=pbc3 OUTEST=esti0;
    MODEL followup*status(0)=tment bililarge logbilipos;
RUN;

DATA esti0; SET esti0;
     IF _NAME_='followup' THEN DO;
        bnoint=tment; key=1;
        OUTPUT;  
     END;
RUN;

DATA esti0; SET esti0 (KEEP=bnoint key); RUN;
    
DATA esti; MERGE esti esti0; BY key; RUN;

DATA plot; SET esti;
           DO k=1 to 5000; i=k/10;
              logi=log(i); 
              linpred=btment+blp*(i>=17.1)+bpp*(i>=17.1)*(logi-log(17.1));
              var=vt+vl*(i>=17.1)+vp*(i>=17.1)*(logi-log(17.1))**2
                  +2*ctl*(i>=17.1)+2*ctp*(i>=17.1)*(logi-log(17.1))
                  +2*clp*(i>=17.1)*(logi-log(17.1));    
              lower=linpred-1.96*sqrt(var);
              upper=linpred+1.96*sqrt(var);              
              logi=log(i);
              nul=bnoint;
              OUTPUT;
            END;
RUN;

PROC GPLOT DATA=plot;
      PLOT linpred*logi=1 lower*logi=2 upper*logi=2 nul*logi=2
                               /FRAME OVERLAY HAXIS=axis1
                               VAXIS=axis2;
          SYMBOL1 V=NONE INTERPOL=SPLINE L=1;
          SYMBOL2 V=NONE INTERPOL=SPLINE L=33;
          axis1 ORDER=0 TO 6 BY 1
                MINOR=NONE LABEL=('log(Bilirubin)');
          axis2 ORDER=-3 TO 2 BY 1
                MINOR=NONE  LABEL=(A=90 R=0 'Linear predictor');
RUN;