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;

data fever;
set fever; 

if death<0 then delete;
/* Definition below makes intercept correspond to mage=30 */
mage30=mage-30;
RUN;

proc genmod data=fever;
    class parity;
    model episodes=parity mage30 / dist=poisson link=log;
    output out=fit p=predicted reschi=residual;
run;

data b0;
set fit;

/* calculation of poisson probabilities for expected values */ 
episode_group=0; expected=poisson(predicted,0); output;
episode_group=1; expected=poisson(predicted,1)-poisson(predicted,0); output;
episode_group=2; expected=poisson(predicted,2)-poisson(predicted,1); output;
episode_group=3; expected=1-poisson(predicted,2); output;
run;

data b1;
set b0;

indicator=(episodes=episode_group);
if episode_group=3 then indicator=(episodes>2);
run;

proc means nway noprint N sum data=b1;
    class episode_group;
    var indicator expected;
output out=chi2 N=antal sum=Sum_Observed sum_Expected;
run;


data chi2;
set chi2;

deviation=(Sum_Observed-sum_Expected)/sqrt(sum_Expected);
deviation_squared=deviation**2;
run;

PROC FORMAT;
    VALUE grouplabels
        0 = "0"
        1 = "1"
        2 = "2"
        3 = "3+";
RUN;

/*************************************************************
    NOTE: The last entry in Table 7.2.4 should be 4.189, not 17.545
    17.545 is the squared deviation
    The chi-square test is correct
    *******************************************************/
    
title 'Note: the last entry in Table 7.2.4 should be 4.189, not 17.545, see NOTE in program';
proc print data=chi2;
 var episode_group Sum_Observed sum_Expected deviation deviation_squared;
 FORMAT episode_group grouplabels. ;
run;

title 'Chi-square test is the sum of the squared deviations, 19.97';
proc means n sum data=chi2;
  var deviation deviation_squared;
run;