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

/*
Creating data set "sums", which is a cross table of death vs alcogroups.
*/

PROC MEANS DATA = fever SUM N NOPRINT;
    CLASS alcogroup;
    VAR death;
    OUTPUT OUT = sums N=total Sum=death;
RUN;

/*
The model is fitted with no intercept to obtain log(odds),
which are then merged onto the original dataset.
*/

PROC GENMOD DATA = sums;
    CLASS alcogroup;
    MODEL death/total = alcogroup /
        DIST = BINOMIAL
        LINK = LOGIT
        NOINT; /* no intercept */
    OUTPUT OUT = estimates XBETA=logodds;
RUN;

PROC SORT DATA = fever;
    BY alcogroup;
RUN;

PROC SORT DATA = estimates;
    BY alcogroup;
RUN;    

DATA fever;
    MERGE fever estimates (KEEP = logodds alcogroup);
    BY alcogroup;
    RUN;

/*
Mean alcohol consumption within alcogroups is added to the dataset.
*/

PROC MEANS DATA = fever NOPRINT MEAN;
    CLASS alcogroup;
    VAR alco;
    OUTPUT OUT = means
        MEAN = alcomean ;
RUN;

DATA means;
    SET means;
    IF alcogroup NE . ;
RUN;

DATA fever;
    MERGE fever means;
    BY alcogroup;
RUN;

/*
Linear model with alcomean as explanatory variable.   
*/

PROC GENMOD DATA = fever DESCENDING;
    MODEL death = alcomean /
        DIST = BINOMIAL
        LINK = LOGIT;
    OUTPUT OUT = est
        XBETA = newlogodds;
RUN;

/*
log(odds) plotted against mean consumption of alcohol (within alcogroups).
Added line corresponding to the fitted model above.
*/

PROC GPLOT UNIFORM DATA = est;
    PLOT logodds*alcomean=1 newlogodds*alcomean=2/ OVERLAY
        HAXIS=AXIS1 VAXIS=AXIS2; 
    AXIS1
        MINOR=NONE
        ORDER=(0 TO 5 BY 1)        
        LABEL=('Drinks per week');
    AXIS2
        MINOR=NONE
        ORDER=(-5.0 TO -4.0 BY .2)
        LABEL=(A=90 R=0 'log(odds)');
    SYMBOL1 VALUE=CIRCLE COLOR=BLACK;
    SYMBOL2 VALUE=NONE INTERPOL=JOIN COLOR=BLACK;
RUN;