/*********************************************************
Number of survivors following the sinking of the titanic
classified by age, economic status, and gender.
*********************************************************/
data grouped;
input age $ sex $ status $ surv n;
adult=(age='Adult');
male=(sex='Male');
second=(status='Second');
third=(status='Third');
crew=(status='Crew');
datalines;
Adult  Male    First   57  175
Adult  Male    Second  14  168
Adult  Male    Third   75  462
Adult  Male    Crew   192  862
Child  Male    First    5    5
Child  Male    Second  11   11
Child  Male    Third   13   48
Child  Male    Crew     0    0
Adult  Female  First  140  144
Adult  Female  Second  80   93
Adult  Female  Third   76  165
Adult  Female  Crew    20   23
Child  Female  First    1    1
Child  Female  Second  13   13
Child  Female  Third   14   31
Child  Female  Crew     0    0
;
run;

/***********************************************************
Now rewrite the dataset so that there is one observation for
every individual. (surv=0 or 1 and n=1 for every observation)
************************************************************/
data individ;
set grouped;
obs=n;
r=surv;
do i=1 to obs;
surv=(i le r);
n=1;
output;
end;
return;
run;

proc logistic data=grouped;
title 'Model fitted to grouped data';
model surv/n = second third crew male
               / influence iplots;
run;

proc logistic data=grouped;
title 'Model fitted to grouped data';
model surv/n = second third crew male
               / influence iplots;
output out=diag difchisq=difchisq difdev=difdev pred=pred
                reschi=reschi resdev=resdev h=h;
run;

/* set the graphics environment */
goptions reset=all gunit=pct ftext=hwpsl009 /* helvetica */
         htitle=5.0 htext=4.0 noprompt
         rotate=portrait horigin=1in vorigin=0.5in
         vsize=12.0cm hsize=12.0cm
         device=pslepsf gsfname=grfout gsfmode=replace
         device=win target=winprtg ftext=swiss
;

proc gplot data=diag;
title 'Change in deviance vs predicted probability';
axis1 label=(angle=-90 rotate=90 'Dev. change');
plot difdev*pred / vaxis=axis1;
run;

title 'Deviance residual vs predicted probability';
axis1 label=(angle=-90 rotate=90 'Dev. residual');
bubble resdev*pred=h  / vaxis=axis1;
run;
quit;

proc logistic data=individ;
title 'Individual data, casewise approach (the default)';
model surv/n = second third crew male
               / scale=none;
run;

proc logistic data=individ;
title 'Individual data, collapsing approach';
model surv/n = second third crew male
               / scale=none aggregate;
run;

proc logistic data=individ;
title 'Individual data, contingency table approach';
model surv/n = second third crew male
               / scale=none aggregate=(second third crew male adult);
run;

proc logistic data=grouped;
title 'Grouped data, contingency table approach (the default)';
model surv/n = second third crew male
               / scale=none;
run;

proc logistic data=individ;
title 'Individual data, Hosmer-Lemeshow test';
model surv/n = second third crew male
               / lackfit;
run;








proc genmod data=grouped;
class status sex;
model surv/n = status sex / dist=bin;
run cancel;