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;

DATA pbc3; SET pbc3;
    bilisq=bili*bili;
    fail=(status>0);
    PROC SORT; BY followup;
RUN;

PROC PHREG DATA=pbc3 OUTEST=esti COVOUT;
    MODEL followup*status(0)=bili bilisq;
    OUTPUT OUT=diagn DFBETA=diag1 diag2/ORDER=DATA;
RUN;

DATA pbc3; MERGE pbc3 diagn; BY followup; key=1; run;

DATA esti; SET esti;
key=1;
IF _NAME_='followup' THEN DELETE;
IF _NAME_='bili' THEN DO; v1=bili; c=bilisq; END;
IF _NAME_='bilisq' THEN DO; v2=bilisq; END;
RETAIN v1 c;    
OUTPUT;
RUN;



DATA esti; SET esti (KEEP=v1 v2 c key);
IF v2=. THEN DELETE;
RUN;
    
DATA pbc3; MERGE PBC3 esti; by key; run;    

/* compute Cook's distance based on deletion diagnostics and parameter variance */
    
DATA pbc3; SET pbc3;
    cook=(diag1*(diag1*v2-c*diag2)+diag2*(-diag1*c+diag2*v1))/(2*(v1*v2-c*c));
RUN;

%begin_multilayout(1,2);

PROC GPLOT DATA=pbc3 GOUT=plots; 
    PLOT cook*bili=fail/ 
        HAXIS=axis1 VAXIS=axis2;
    axis1 ORDER=0 TO 500 BY 100 MINOR=NONE
        LABEL=('Bilirubin');
    axis2 ORDER=0 TO 0.25 BY 0.05
        MINOR=NONE LABEL=(A=90 R=0 'Cook');
    symbol1 V=circle;
    symbol2 V=PLUS;
RUN;

PROC GPLOT DATA=pbc3  GOUT=plots; 
    PLOT cook*followup=fail/
        HAXIS=axis1 VAXIS=axis2;
    axis1 ORDER=0 TO 6 BY 1 MINOR=NONE
        LABEL=('Years');
    axis2 ORDER=0 TO 0.25 BY 0.05
        MINOR=NONE LABEL=(A=90 R=0 'Cook');
    symbol1 V=circle;
    symbol2 V=PLUS;
RUN;

%end_multilayout(1,2);