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;

/* Compute fractional polynomials for powers q=-0.5 and -2
    plus interactions with tment*/

DATA pbc3; SET pbc3;
    bili_05=1/sqrt(bili); bili_05pred=bili_05*tment;
    bili_2=1/bili**2; bili_2pred=bili_2*tment;
RUN;

PROC PHREG DATA=pbc3;
    MODEL followup*status(0)=tment bili_05 bili_2 bili_05*tment bili_2*tment;
    ODS OUTPUT PARAMETERESTIMATES=estimates(KEEP=ESTIMATE);
RUN;

PROC TRANSPOSE DATA=estimates OUT=esttrans; RUN;

DATA esttrans; SET esttrans
     (RENAME=(col1=btment col2=b_05 col3=b_2 col4=b_05p col5=b_2p));
RUN;

DATA plot; SET esttrans;
           DO k=23 to 4530; i=k/10;
              linpred0=b_2/(i*i)+b_05/sqrt(i);
              linpred1=linpred0+btment+b_2p/(i*i)+b_05p/sqrt(i); ;
              logi=log(i);
              OUTPUT;
              END;
RUN;

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