%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter6;
PROC IMPORT OUT= WORK.dwaine 
            DATAFILE= "&location\chapter6.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH06FI05$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;

proc print data=dwaine;
title 'listing of data refer to figure 6.5 p. 237';run;
proc reg data=dwaine;
model sales=children income/stb;
title 'regression analysis of data refer to figure 6.5 p. 237';run;
proc glm data=dwaine;
model sales=children income/inverse xpx;
title 'analysis of variance summary table, x''x, inv(x''x). refer figure 6.5 p. 237';run;
proc g3d data=dwaine;
scatter income*children=sales/grid noneedle shape='balloon';
title '3d scatterplot refer to Figure 6.6';run;
proc g3d data=dwaine;
scatter income*children=sales/grid noneedle rotate=0 shape='balloon';
title 'unrotated 3d scatterplot refer to Figure 6.6';run;
proc g3d data=dwaine;
scatter income*children=sales/grid noneedle rotate=30 shape='balloon';
title 'rotated 3d scatterplot refer to Figure 6.6(a)';run;
proc g3d data=dwaine;
scatter income*children=sales/grid noneedle rotate=-60 shape='balloon';
title 'rotated 3d scatterplot refer to Figure 6.6(b)';run;
proc g3d data=dwaine;
scatter income*children=sales/grid noneedle rotate=0 to 360 by 30 shape='balloon';
title 'multiple rotations refer to Figure 6.6 Exploration';run;
proc means data=dwaine;
title 'summary statistics for data set- needed to make plot';run;
/*Make a variable which contains the observation number*/
data dwaine;set dwaine;
id+1;
/*Make a data set of possible values for the predicted regression surface*/
data grid;
sales=.;
do children=38 to 92 by 1;
 do income=15 to 20 by .2;
 output;
 end;
 end;
data dwaine2;set dwaine grid;
proc reg data=dwaine2;
model sales=children income/stb;
output out=results
p=yhat r=residual rstudent=rstudent h=hatvalue;
id id;
data resultsp;set results;if sales=.; sales=yhat;
 proc g3d data=resultsp;
 plot income*children=sales/rotate=0 to 360 by 30 ctop=red cbottom=green
xticknum=3 yticknum=3 ;
title 'predicted response surface refer to figure 6.7';run;

data fit;set results end=last;if sales=.; sales=yhat;

if last;  * Create two 'dummy' observations with the min and max x, y and z values ;
  children=38;
  income=15;
  sales=135;
colorval='white'; shapeval='point';  *  'Hide' the 'dummy' points ;
  output;
children=92;
  income=20;
  sales=250;
colorval='white'; shapeval='point';  *  'Hide' the 'dummy' points ;
  output;
run;
data needles;set results;if sales^=.;
length caseid text $12;
caseid=' ';
if abs(rstudent)>tinv(.9,21-3,0)
 |hatvalue>2*3/21 then caseid=put(id,3.);
 xsys='2';ysys='2';zsys='2';
 y=income;x=children;
 z=yhat;function='MOVE    ';output;
 if residual>0 then 
   do;
   line=1;color='BLACK';position='2';
   end;
 else 
 do;
 line=1;color='RED';position='8';
 end;
z=sales;function='DRAW    ';output;
function='SYMBOL';text='+';output;
function='LABEL';text=caseid;output;

 proc g3d data=fit;
 plot income*children=sales/rotate=0 to 360 by 90 ctop=red cbottom=green
   xticknum=3 yticknum=3 anno=needles grid;
title 'plot of estimated regression surface refere to figure 6.7';run;

data fitresult;set fit resultsp;
 proc g3d data=fitresult;
 plot income*children=sales /rotate=0 to 360 by 60 ctop=yellow cbottom=blue
   xticknum=3 yticknum=3 anno=needles grid;
title 'plot of estimated regression surface refere to figure 6.7';run;
proc reg data=dwaine;
model sales=children income/stb;
plot residual.*predicted. residual.*children residual.*income;
output out=checkit
p=yhat r=residual rstudent=rstudent h=hatvalue;
id id;
title 'diagnostic plots refer Figure 6.8 p. 242';
data checkit;set checkit;
inter=children*income;
absresid=abs(residual);
symbol i=rlcli95 v=dot;
proc gplot;
plot residual*inter;
title 'diagnostic plot of interaction. Refer Figure 6.8 p. 242';run;
proc gplot;
plot absresid*yhat;
title 'additional diagnostic plots. Refer Figure 6.8(a) p. 243';run;
proc capability data=checkit;var residual;
qqplot residual;
title 'normal probability plot';run;
data checknorm;set checkit;
normr=residual;
proc rank normal=blom data=checknorm out=checknorm;var normr;
proc corr;var normr residual;
title 'calculation of correlation between observed and normal distribution errors refer p. 243';run;
proc glm;
model sales=children income/solution;
title 'anova summary statistics refer to page 244.';run;
proc reg;
model sales=children income/corrb covb alpha=.05 clb;
title 'anova summary statistics refer to page 245.';run;
data newdata;
children=65.4;income=17.6;output;
data dwaine2;set newdata dwaine;
proc reg data=dwaine2;
model sales=children income/alpha=.05 clm;
title 'anova summary statistics refer to page 245.';run;
data newdata;
children=65.4;income=17.6;output;
children=53.1;income=17.7;output;
title 'looking at prediction limits for new observations';
data dwaine2;set newdata dwaine;
proc reg data=dwaine2;
model sales=children income/alpha=.05 cli;
title 'anova summary statistics refer to page 245. Bonferroni correction';run;

proc glm data=dwaine2 outstat=stats;
model sales=children income;
output out=scheffeb p=psales stdi=stdi;
title 'calculating Scheffe interval Refer page 247';
run;
data stats;set stats;
if _source_='ERROR';
drop _name_ _type_ f prob;
grp=1;
data scheffeb;set scheffeb;grp=1;
data scheffeb;merge scheffeb stats;by grp;
drop _source_ grp SS;
if sales=.;
*What are our conditions?;
alpha=.10;
predicts=2;
Label
smult= 'Scheffe Multiple'
bmult= 'Bonferroni Multiple'
slow= 'Scheffe Lower Bound'
shigh= 'Scheffe Upper Bound'
blow= 'Bonferroni Lower Bound'
bhigh= 'Bonferroni Upper Bound';

smult=sqrt(2*finv(1-alpha,2,df));
bmult=tinv(1-alpha/(2*predicts),df);
slow=psales-smult*stdi;
shigh=psales+smult*stdi;
blow=psales-bmult*stdi;
bhigh=psales+bmult*stdi;
proc print data=scheffeb;run;
data _null_;set scheffeb;
file print;
put 'Assuming a family confidence coefficient of 1-' alpha 4.2/
    @5 'The Lower and Upper Scheffe intervals are: ' slow 8.2 ' and 'shigh 8.2 ' respectively, and'/
	@5 'The Lower and Upper Bonferroni intervals are: ' blow 8.2 ' and ' bhigh 8.2 ' respectively'/
	@10 '(Assuming a Scheffe multiple of ' smult 12.8 ' and a Bonferroni multiple of ' bmult 12.8 ')';
	put 'assuming error df of ' df;
title 'Scheffe and Bonferroni Bounds on multiple new observations. Refer to Page 247';run;