%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter20;
PROC IMPORT OUT= WORK.premium 
            DATAFILE= "&location\chapter20.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH20TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value regfmt
 1="1=East"
 2="2=West";
value sizefmt
 1="1=Small"
 2="2=Medium"
 3="3=Large";
data premium;set premium;
format region regfmt. size sizefmt.;
proc means;
title "Insurance premium example. Refer to Table 20.2a, page 862";run;
axis1 minor=none offset=(5) order=100 to 220 by 20;
axis2 minor=none offset=(5);
symbol1 v=dot i=join;
proc gplot;
plot premium*size=region/noframe haxis=axis2 vaxis=axis1;
title "Alternative to Figure 20.2, page 882";run;
axis1 minor=none offset=(5) order=100 to 220 by 20;
axis2 minor=none offset=(5) order=1 2;

proc gplot;
plot premium*region=size/noframe haxis=axis2 vaxis=axis1 ;
title "Refer to Figure 20.2, page 882";run;

proc glm;class size region;
model premium=size|region;
title "why you can't do glm with N=1";run;
proc glm data=premium;class size region;
model premium=size region;
output out=outprem p=predprem;
title "using the interaction of size with region as error term";
title2 "Refer to page 883-884";run;
proc print data=outprem;var size region predprem;
title "Estimated Treatment Means. refer to pp. 884-885";run;
proc glm;class size region;
model premium=size|region;
test h=size region e=size*region; 
title 'same test but using test statement';run;
proc means data=premium;class region size;var premium;
title "Observed Cells- Refer to pp. 884-885";run;

ods listing close;
proc glm data=premium;
  class region size;
  model premium = region size;
  ods output overallanova=overall modelanova=model;
run;
/*Tukey test for additivity*/
/*Here are two ways to do this. The first captures the statistics described in the text and then
  calculates the Tukey test. The second (and much easier) approach invokes the twoway.sas macro written
  by Michael Friendly and, in addition to calculating the test, also produces useful graphics for interpreting
  whether the nonadditivity could be due to other properties of the data set*/
quit;
ods listing;
ods output close;
data _null_;
  set overall;
  if source='Corrected Total' then call symput('overall', ss);
run;
data _null_;
  set model ;
  if hypothesistype=1 and source='size' then call symput('ssa', ss);
  if hypothesistype=1 and source='region' then call symput('ssb', ss);
  if hypothesistype=1 and source='size' then call symput('dfa', df);
  if hypothesistype=1 and source='region' then call symput('dfb', df);
run;
proc means data=premium;var premium;
output out=meant mean=meant;
title 'calculate grand mean';run;
data _null;set meant;
call symput('meant',meant);run;
proc means data=premium;class region;var premium;
output out=meanr mean=meanr;
title 'calculate means by region';run;
proc means data=premium;class size;var premium;
output out=means mean=means;
title 'calculate means by size';run;
proc sort data=means out=means;by size;
proc sort data=premium out=premium;by size;
data allmeans;merge premium means;by size;
proc sort data=meanr out=meanr;by region;
proc sort data=allmeans out=allmeans;by region;
data allmeant;merge allmeans meanr;by region;
meant=&meant;
if region>. and size>.;
drop _type_ _freq_;
if region>. and size>.;
ss=(meanr - meant)*(means - meant)*premium;
proc means data=allmeant;var ss;
output out=totss sum=sumss;
title 'calculate sstotal';
data _null_;set totss;
call symput('total',sumss);run;

data final;
  msa = &ssa/(&dfb+1);
  msb = &ssb/(&dfa+1);
  ssab = (&total*&total) / ( msa*msb );
  ssrem = &overall - &ssa - &ssb - ssab;
  f = ssab/( ssrem/((&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1)) );
  p_value = 1- cdf('F',f,  1, (&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1) );
  adddendf=(&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1);
  addnumdf=1;
run;
proc print data=final;
title 'Tukey additivity test';run;
run;
options macrogen;
%include "&location\twoway.sas";

%twoway(data=premium,var=size region,response=premium);