%include "C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter8\3dplotscatter.sas";
%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter8;
data fig8_1;
do x=-3 to 3 by .1;
model="a";
y=52+(8*x)-(2*x*x);
output;
model="b";
y=18-(8*x)+(2*x*x);
output;
end;
proc sort data=fig8_1;by model;
symbol1 c=green v=none l=1 i=Join;
symbol2 c=red v=none l=10 i=Join;
proc sort data=Fig8_1 out=Fig8_1;by model;
proc gplot;
plot y*x/vaxis=0 to 60 by 10 
             haxis=-3 to 3 ;by model;
title "Data from Figure 8.1, page 295";run;
data fig8_2;
do x=-3 to 3 by .1;
model="a";
y=22.45+1.45*x+.15*x*x+.35*x*x*x;
output;
model="b";
y=16.3-1.45*x-.15*x*x-.35*x*x*x;
output;
end;
proc sort data=fig8_2;by model;
symbol1 c=green v=none l=1 i=Join;
symbol2 c=red v=none l=10 i=Join;
proc sort data=Fig8_2 out=Fig8_2;by model;
proc gplot;
plot y*x/vaxis=0 to 30 by 10 
             haxis=-3 to 3 ;by model;
title "Data from Figure 8.2, page 297";run;
data fig8_3;
do x1=-10 to 10 by .1;
 do x2=-10 to 10 by .1;
 y=1740-4*x1*x1-3*x2*x2-3*x1*x2;
 output;
 end;
end;
proc g3d data=fig8_3;
 plot x1*x2=y/rotate=0 to 360 by 30 ctop=red cbottom=green
xticknum=3 yticknum=3 ;
title "Figure 8.3a, p. 298";run;
proc gcontour;
 plot x1*x2=y;
 title "Figure 8.3b, p. 298";run;

PROC IMPORT OUT= WORK.power 
            DATAFILE= "&location\chapter8.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH08TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Refer to Table 8.1, page 300)";run;
data power;set power;
id+1;
label 
chargec="centered version of charge"
tempc="centered version of temp";
chargec=charge;
tempc=temp;
proc standard data=power out=power mean=0;var chargec tempc;
data power;set power;
chargec=chargec/.4;
tempc=tempc/10;
proc print;
title "Refer to Table 8.1, page 300";run;
data power;set power;
/* centered product and quadratic terms*/
label chargec2="quad. centered charge"
      tempc2="quad. centered temp"
	  chc_tpc="int. centered charge x temp";
chargec2=chargec*chargec;
tempc2=tempc*tempc;
chc_tpc=chargec*tempc;
/* raw product and quadratic terms*/
label charge2="quad. raw charge"
      temp2="quad. raw temp"
	  ch_tp="temp x charge";
charge2=charge*charge;
temp2=temp*temp;
ch_tp=charge*temp;
proc print;
title "Refer to Table 8.1, page 300";run;
proc corr;var temp temp2 charge charge2 ch_tp;
title "correlations between raw predictors. Refer to Page 301";run;
proc corr;var tempc tempc2 chargec chargec2 chc_tpc;
title "correlations between centered predictors. Refer to Page 301";run;
proc reg;
model cycles= temp temp2 charge charge2 ch_tp;
title "regression analysis based on raw variables";run;
proc reg;
model cycles=tempc tempc2 chargec chargec2 chc_tpc/stb;
output out=powerc r=rcycles;
plot residual.*predicted. residual.*chargec residual.*tempc nqq.*predicted.;
title "regressions results. Refer Figure 8.4, p. 302, Figure 8.5, page 303";run;
proc capability data=powerc;var rcycles;
qqplot rcycles;
title "normal probability plot - Reference Figure 8.4, p. 302";run;
proc glm data=power;
model cycles=tempc tempc2 chargec chargec2 chc_tpc/solution;
title "reduced model: Refer p. 301-302 & Figure 8.4";run;

proc glm;class chargec tempc;
model cycles=chargec|tempc;
title "full model: Refer p. 302";run;
data _null_;
ndf=3;
ddf=2;
f_lackf=((5240.43860-1404.66667)/3)/(1404.66667/2);
critf=finv(1-.05,ndf,ddf);
probf=probf(f_lackf,ndf,ddf);
file print;
put "F test for lack of fit " f_lackf 8.4 /
 " has probability " probf 8.4 /
 " with critical F of " critf 8.4 " with " ndf " and " ddf "degrees of freedom";
 run;

/*The following code shows the differences in coding in the solution with and without the class statement*/
proc glm;
model cycles=tempc tempc2 chargec chargec2 chc_tpc/solution;
title "reduced model: Refer p. 301-302 & Figure 8.4";run;

proc glm;class chargec tempc;
model cycles=chargec|tempc/solution;
title "illustration of solution option with class statement";run;

proc reg;
model cycles=tempc tempc2 chargec chargec2 chc_tpc/stb;
quadint: test chargec2=tempc2=chc_tpc=0;
title "Test of whether quadratic and interaction terms are all zero Refer to pp. 303-304";run;
proc reg;
model  cycles=charge temp /stb;
title "Raw score main effect model. Refer to Equation 8.19, p. 304";run;
/*Here"s a macro to do those 3d plots*/
/*%macro scatter3d(data, /*Name of data set*/
/*                 x, /*Name of X variable (1st predictor)*/
/*                 y, /*Name of Y variable (2nd predictor)*/
/*                 z, /*name of criterion variable*/
/*              minx, /*Minimum X value (1st predictor)*/
/*              miny, /*Minimum Y value (2nd predictor)*/
/*              minz,  /*Minimum z (criterion) value*/
/*              maxx, /*Maximum X value (1st predictor)*/
/*              maxy, /*Maximum Y value (2nd predictor)*/
/*              maxz, /*maximum z (criterion) value*/
/*                 n, /*sample size*/
/*              stepx=1, /*grid steps of x*/
/*              stepy=1); /*grid steps of y*/  
data power2;set power;
%scatter3d(data=power2, x=charge, y=temp, z=cycles, minx=.7, miny=12, minz=0, maxx=1.3, maxy=28, maxz=400, n=11, stepx=.1, stepy=1); 
proc reg data=power;
model  cycles=charge temp /alpha=.05 clb;
title "Raw score main effect model. Refer to p. 305";run;
proc iml;
use power;
read all into x;
print x;
preds=x[,{2 3}];
print preds;
pred1=x[,{2}];
pred1t=t(pred1);
orpol1=orpol(pred1,2);
print orpol1;
pred2=x[,{3}];
pred2t=t(pred2);
orpol2=orpol(pred2,2);
print orpol2;
orpreds=orpol1[,{2 3}]||orpol2[,{2,3}];
newp=x[,{1}]||orpreds[,];
print newp;
names={"cycle" "ch1" "ch2" "temp1" "temp2"};
create poworth from newp[colname=names];
append from newp;
proc print data=poworth;
title "orthogonal polynomials for Power Cells- Refer to page 305";run;
proc reg;
model cycle=ch1 ch2 temp1 temp2/stb;
proc corr;
title "note that orthogonal polynomial coefficients are uncorrelated.";
title2 "Refer to discussion on page 305";run;
data Fig_8_7;
do x1=1 to 10 by .01;
figure="8.7a";
x2=3;y=25+2*x1;output;
x2=1;y=15+2*x1;output;
figure="8.7b";
x2=3;y=25+3.5*x1;output;
x2=1;y=15+2.5*x1;output;
figure="8.7c";
x2=3;y=25+.5*x1;output;
x2=1;y=15+1.5*x1;output;
end;
symbol1 c=green v=none l=1 i=rl;
symbol2 c=red v=none l=10 i=rl;
proc sort data=Fig_8_7 out=Fig_8_7;by figure;
proc gplot;
plot y*x1=x2/vaxis=0 to 60 by 15 
             haxis=0 to 10 by 5;by figure;
title "Additive, Reinforcement, and Interference examples. Refer p. 307";run;
data fig8_8;
do x1=0 to 50 by 1;
 do x2=0 to 50 by 1;
 y=10+2*x1+5*x2;
 model="8.8a";
 output;
 y=10+2*x1+5*x2+5*x1*x2;
 model="8.8b";
 output;
 model="8.8c";
 y=10+2*x2+5*x2-5*x1*x2;
 output;
 end;
end;
proc sort data=fig8_8;by model;
proc g3d data=fig8_8;by model;
 plot x1*x2=y/rotate=0 to 360 by 30 ctop=red cbottom=green
xticknum=3 yticknum=3 ;
title "predicted response surface refer to figures 8.8a-8.8c p. 310";run;
proc gcontour;by model;
 plot x1*x2=y;
 title "Contour Plots refer to Figures 8.8d-8.8f, p. 310";run;
data fig8_9;
do x1=-1 to 1 by .1;
 do x2=-1 to 1 by .1;
 y=65+3*x1+4*x2-10*x1*x1-15*x2*x2+35*x1*x2;
 output;
 end;
end;
proc g3d data=fig8_9;
 plot x1*x2=y/rotate=0 to 360 by 60 ctop=red cbottom=green
xticknum=3 yticknum=3 ;
title "predicted response surface Quick Bread Example refer to figure 8.9a p. 311";run;
proc gcontour;
 plot x1*x2=y;
 title "Contour Plots Quick Bread Example refer to Figure 8.9b, p. 311";run;
data fig8_10;set fig8_9;
if x1>.99 or x1=-1;
/*Note use of line types in the symbol statement. Nice values of lines are 1, 2, 26, & 20
or 1, 3, 4 & 8
 (in that order) See the l= option under the symbol statement for the possibilities
 if you need more options */
symbol1 i=rq v=none c=red l=1;
symbol2 i=rq v=none c=blue l=2;
proc gplot;
plot y*x2=x1;
title "Conditional Effect Plot Refer to Figure 8.10, p. 311";run;
PROC IMPORT OUT= WORK.bodyfat 
            DATAFILE= "&location\chapter7.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH07TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print data=bodyfat;
title "Hi, remember me? I'm the Bodyfat Data  refer to Table 7.1";run;
data bodyfat;set bodyfat;
tri_thigh=triceps*thigh;
tri_mid=triceps*midarm;
thigh_mid=thigh*midarm;
proc corr;var triceps thigh midarm tri_thigh tri_mid thigh_mid;
title "correlations of predictors and interaction terms if not centered";run;
proc standard data=bodyfat out=bodyfatc mean=0;var triceps thigh midarm;
data bodyfatc;set bodyfatc;
tri_thigh=triceps*thigh;
tri_mid=triceps*midarm;
thigh_mid=thigh*midarm;
proc corr;var triceps thigh midarm tri_thigh tri_mid thigh_mid;
title "correlations of predictors and interaction terms if centered";run;
proc glm data=bodyfatc;
model bodyfat=triceps thigh midarm tri_thigh tri_mid thigh_mid;
title "extra sums of squares for centered bodyfat data";run;
proc glm data=bodyfat;
model bodyfat=triceps thigh midarm tri_thigh tri_mid thigh_mid;
title "extra sums of squares for uncentered bodyfat data had we not centered";
title2 "Note that Type I SS are same, but, as expected, Type III SS are different";run;
proc reg data=bodyfatc;
model bodyfat=triceps thigh midarm tri_thigh tri_mid thigh_mid;
noint: test tri_thigh=tri_mid=thigh_mid=0;
title "test of no interaction effects. Refer to page 313";run;

data critf;
file print;
critf=finv(1-.05,3,13);put "Critical F value  at alpha=.05 is " critf= " with 3 and 13 df";
title "'Required critical F for test of interaction effects....' refer page 313";run;
PROC IMPORT OUT= WORK.insure 
            DATAFILE= "&location\chapter8.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH08TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class type;
title "Insurance Innovation data- refer to Table 8.2, page 317";run;
/*Labeling your formats can make your graphs easier to understand and can document your
  analyses*/
proc format;
value typef
 0="0=Mutual"
 1="1=Stock";
data insure;set insure;
format type typef.;
proc reg;
model months=size type/stb;
title "regression predicting number of months elapsed from size of firm and type. Refer to page 317, Table 8.3";run;
symbol1 i=rl v=circle c=red;
symbol2 i=rl v=dot c=blue;
proc gplot;
plot months*size=type;
title "plot of regression functions by type of firm, Refer to Figure 8.12, page 318";run;
proc reg;
model months=size/stb;
output out=insure r=resid p=pred;
proc gplot;
plot resid*pred=type;
title "diagnostic plots of informative residuals. This foreshadows Figure 8.16, page 331";run;
data insure;set insure;
size_typ=size*type;
proc reg;
model months=size type size_typ;
title "test of inequality of slopes, refer to page 326 & Table 8.4, p. 327";run;
proc glm;class type;
model months=size|type;
title "same analyses, but done using proc glm refer to page 326 and Table 8.4, p. 327";run;
data _null_;
file print;
critt=tinv(1-.05/2,16);
put "Critical T value for two-tailed T-test at alpha=.05 " critt=;
title "critical T value for test of interaction effect. Refer to page 327";run;
PROC IMPORT OUT= WORK.soap 
            DATAFILE= "C:\Documents and Settings\woodph\My Documents\glm
wood\kutnersolutions\Chapter8\chapter8.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH08TA05$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class line;
title "Summary statistics- Soap Production Referto Table 8.5, p. 330";run;
symbol1 i=rl v=dot c=red;
symbol2 i=rl v=circle c=blue;
proc format;
value linef
 0="0=Line 2"
 1="1=Line 1";

data soap;set soap;
format line linef.;
proc gplot;
plot scrap*speed=line;
title "plot of speed by scrap as a function of Production line, Refer to Figure 8.16, p. 331";run;
proc glm;class line;
model scrap=speed|line/solution;
title "regression solution via glm= Refer to Page 331. Note reflections due to scaling of line";run;
data soap;set soap;
speedlin=speed*line;
proc reg;
model scrap=speed line speedlin/stb corrb covb;
output out=soap r=rscrap p=pscrap;
title "regression solution via reg. Refer to Page 331 & Table 8.6a, p. 332";run;
proc gplot;
plot rscrap*pscrap=line;
title "plots of residual by predicted scrap as a function of line= refer to Fogire 8/17, p. 332";run;
proc sort data=soap out=soap;by line;

/*Brown Forsythe Test- Since speed is not a class variable, the BF test must be done by 
accumulating the relevant statistics and plugging the values in*/
proc reg data = soap ;by line;
  model scrap = speed;
  output out=rsoap r=resid;
title "calculate and save residuals. Refer to page 332";
run;
proc means data = rsoap noprint;
class line;
  var resid;
  output out=medout median=medianr;
title "calculate median of residuals within class";
run;
proc sort data=rsoap out=rsoap;by line;
proc sort data=medout out=medout;by line;
title "merge residuals with medians & calculate d. Refer to page 332";
data mtemp;
  merge rsoap medout;
  by line;
  d = abs(resid - medianr);
title "calculate absolute value of residual and median. Refer to page 332";
run;

proc ttest data = mtemp;
  class line;
  var d;
title "Brown-Forsythe results. Refer to page 333";
run;

proc reg data=soap;
model scrap=speed line speedlin;
mainef: test line=speedlin=0;
inter: test speedlin=0;
title "test of interaction between line and speed as well as between two production lines with interaction. Refer to Page 333";run;
proc glm data=soap;
model scrap=speed|line;
contrast "int" speed*line 1 ;
contrast "line&int" speed*line 1
                    line 1 ;
title "test of differences in slope using Proc GLM with contrast statements. Refer to Page 333";run;
/*Generating some data for exxample on page 334- "Pay no attention to that man behind the curtain*/
data startvals;
do i=1 to 50;
x=ranuni(919821);
error=rannor(28391);
machine=0;
output;
machine=1;
output;
end;
run;
data calibrat;set calibrat;drop error;*Need to do this if we're going to have residual have the same name;
proc reg noprint;
model error=x;
output out=calibrat r=error;
proc standard data=calibrat out=calibrat mean=0 std=1;var error;
proc standard data=calibrat out=calibrat mean=0 std=1;var x;
data calibrat;set calibrat;
y=.1+.2*x+.3*x*x+.4*machine-.5*machine*x-.6*machine*x*x+error;
x2=x*x;
xm=x*machine;
x2m=x*x*machine;
proc reg;
model y=machine x x2 xm x2m;
mach: test machine=xm=x2m=0;run;
proc glm;
model y=x x*x machine machine*x machine*x*x/solution;
title "simulated data for the analysis presented on page 334";run;