%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter17;
PROC IMPORT OUT= WORK.kenton 
            DATAFILE= "&location\Chapter17.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH16TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value designf
 1="1=3 colors, with cartoons"
 2="2=3 colors, without cartoons"
 3="3=5 colors, with cartoons"
 4="4=5 colors, without cartoons";
data kenton;set kenton;
format design designf.;
label cases="Cases Sold";
proc means n sum mean;class design;var cases;
title "kenton foods example- means by design Refer Table 17.1, page 734";run;
proc means n sum mean;var cases;
title "kenton foods example- overall means. Refer Table 17.1, page 734";run;
proc glm;class design;
model cases=design;
title "GLM summary table. Refer Table 17.1, page 734";run;
PROC IMPORT OUT= WORK.rust 
            DATAFILE= "&location\Chapter17.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH17TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value brandf
 1="1=A"
 2="2=B"
 3="3=C"
 4="4=D";run;
proc print;
title "Rust inhibitor data. Refer to Table 17.2a, page 735";run;
data rust;set rust;
format brand brandf.;run;
proc glm;class brand;
model inhibit=brand;
title "Rust inhibitor data. Anova summary table. Refer to Table 17.2b, page 735";run;
proc univariate data=kenton;class design;var cases;
histogram cases/normal;
title "means and observed distributions by brand";run;
 
proc gchart data=kenton;
vbar design/type=mean sumvar=cases discrete noframe ref=18.63158 ;
title "Bar Chart of raw means Kenton Foods. Refer to Figure 17.2a, page 736";
title2 "Note that the reference line mentioned as desirable in the text is included in the figure provided here";run;
axis1 order=0 to 35 by 10;
symbol1 i=join c=green;
proc gbarline data=kenton;
bar design/sumvar=cases noframe type=mean ref=18.63158 clm=95 discrete axis=axis1 cerror=blue;
plot /sumvar=cases type=mean raxis=axis1;
title "Bar Charts with Cell Mean Confidence Bars";run;

proc means data=kenton;var cases;
title "finding grand mean";run;
proc means data=kenton;class design;var cases;
output out=kentonm mean=mcases;
symbol1 i=join v=dot;
axis1 order=0 to 30 by 10 ;
axis2 offset=(3 cm);
proc gplot data=kentonm;
plot mcases*design/vaxis=axis1 haxis=axis2 vref=18.63 cvref=red noframe;
title "Main Effects Plot of raw means of Kenton Foods. Refer to Figure 17.2b page 736";run;
axis1 order=5 to 35 by 10 ;
proc gchart data=kenton;
vbar design/clm=95 type=mean sumvar=cases discrete axis=axis1 noframe;
title "Bar Chart of Raw Means Kenton Foods showing individual cell confidence intervlas";
footnote1 "this is *not* equal to Figure 17.3 because individual cells standard deviations are used for";
footnote2 "computation of the confidence intervals but such charts of sometimes useful for comparison with";
footnote3 "confidence intervals from GLM's outlined below";run;
/*Figure 17.3, page 739: Ods is used to capture the lsmeans from glm. 
  These can be quickly plotted with timeplot for typewriter graphics or, 
  more aesthetically, can be plotted with gplot. (I'm unsure how to generate the bar-interval
  graph of figure 17.3a on page 739 in SAS. Your easiest option here is ot export the
  data set of LSmeans (e.g., data set temp below) to excel and make the bar chart with error bar from there.*/
footnote1 ;footnote2 ;footnote3 ;
proc glm data=kenton;class design;
model cases=design;
lsmeans design/cl;
ods output LSMeancl=temp;
run;quit;
proc print data=temp;
title "captured LS means from GLM";run;

proc timeplot;
plot (lowercl lsmean uppercl)=design/overlay hiloc;
title "Interval plots of captured LS means and CI's via timeplot. Refer to Figure 17.3b, page 739";run;

symbol1 interpol=hiloctj
       cv=blue
       ci=red
       width=4;

axis2 offset=(3 cm);
data hiloval;set temp;
value=lowercl;output;
value=lsmean;output;
value=uppercl;output;
label value="Cases Sold";

proc gplot data=hiloval;
plot value*design/overlay haxis=axis2;
title "Interval plot of LS means and confidence intervals. Refer to Figure 17.3b, page 739";run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
estimate 'd3 vs. d4' design 0 0 1 -1;
title "Test for differences between factor means. Refer to example on Page 740.";
title2 "note that default alpha is .05";
run;
proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
contrast 'd3 vs. d4' design 0 0 1 -1;
title "Equivalent test for differences between factor means using Contrast statement. Refer to example on Page 740.";
title2 "note that no confidence interval for difference is calculated";
run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
estimate '3-color vs. 5-color' design .5 .5 -.5 -.5;
title "test of contrast in factor level means. Refer to example on page 743";run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
estimate 'design 1 vs. all others' design .75 -.25 -.25 -.25;
title "test of first design versus all remaining";run;


proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl;
estimate 'grand mean' intercept 1 design .35 .28 .12 .25;
title "estimate of a linear combination, in this case, the overall mean sales";
title2 "refer to page 744";run;
proc glm data=rust;class brand;
model inhibit=brand;
lsmeans brand/cl alpha=.05 adjust=tukey pdiff tdiff;
ods output LSMeancl=lsrust;
title "Tukey pairwise comparisons simultaneous intervals. Refer to page 748 & Table 17.3, page 749";run;quit;
symbol1 interpol=hiloctj
       cv=blue
       ci=red
       width=4;

axis2 offset=(3 cm);
data hiloval;set lsrust;
value=lowercl;output;
value=lsmean;output;
value=uppercl;output;
label value="Performance Score";

proc gplot data=hiloval;
plot value*brand/overlay haxis=axis2;
title "Interval plot of LS means and confidence intervals. Refer to Figure 17.4, page 748";run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=tukey pdiff tdiff;
ods output lsmeancl=lsfood;
symbol1 interpol=hiloctj
       cv=blue
       ci=red
       width=4;
axis2 offset=(3 cm);

data hiloval;set lsfood;
value=lowercl;output;
value=lsmean;output;
value=uppercl;output;
label value="Cases Sold";
proc gplot data=hiloval;
plot value*design/overlay haxis=axis2;
title "Interval plot of LS means and confidence intervals. Refer to pp. 750-751";run;
proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=scheffe pdiff tdiff;
title "Pairwise comparisons using Scheffe";run;
proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=scheffe pdiff tdiff;
estimate '3 color vs 5 color' design .5 .5 -.5 -.5;
estimate 'w cartons vs w/o cartons' design .5 -.5 .5 -.5;
estimate '3 color carton vs 3 color no carton' design 1 -1 0 0;
estimate '5 color carton vs 5 color no carton' design 0 0 1 -1;
ods output estimates=contrasts overallanova=anovatable;
title "Planned Scheffe contrasts. Refer to pp., 754-755";run;quit;
data _null_;set anovatable;
if source='Model' then call symput('modeldf',df);
if source='Error' then call symput('errordf',df);
data contrasts;set contrasts;
file print;
s2=&modeldf*finv(.90,&modeldf,&errordf);
s=sqrt(s2);
lowerci=estimate-stderr*s;
upperci=estimate+stderr*s;
put 'contrast ' parameter 'Estimate ' estimate best10. 'Std. Error ' stderr best10. 
/ @5 'S-Squared ' s2 10.2 'S ' s 10.2 'Lower CI ' lowerci 10.2 'Upper CI ' upperci 10.2;
run;
proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=sidak pdiff tdiff ;
title "Pairwise comparisons using Sidak";run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=simulate(report) pdiff tdiff ;
title "Pairwise comparisons using Simulate";run;

proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
lsmeans design/cl alpha=.1 adjust=scheffe pdiff tdiff;
estimate '3 color vs 5 color' design .5 .5 -.5 -.5;
estimate 'w cartons vs w/o cartons' design .5 -.5 .5 -.5;
ods output estimates=bcontrasts overallanova=anovatable;
title "Planned Bonferroni contrasts. Refer to pp., 756-757";run;quit;
data _null_;set anovatable;
if source='Model' then call symput('modeldf',df);
if source='Error' then call symput('errordf',df);

data bcontrasts;set bcontrasts;
file print;
conf=.975;
alpha=1-conf;
bonf=tinv(1-alpha/(2*2),&errordf);*Refer to 17.46a, page 756 denominator is 2 x # contrasts considered;
lowerci=estimate-stderr*bonf;
upperci=estimate+stderr*bonf;
put 'Parameter ' parameter ' estimate ' estimate best8. 'Std. Error ' stderr best8. +2 'B ' bonf best8. +2 'lowerci ' lowerci best8. +2 'upperci ' upperci best8.;run;
data scontrasts;set bcontrasts;
file print;
s2=&modeldf*finv(.975,&modeldf,&errordf);
s=sqrt(s2);
lowerci=estimate-stderr*s;
upperci=estimate+stderr*s;
put 'contrast ' parameter 'Estimate ' estimate best10. 'Std. Error ' stderr best10. 
/ @5 'S-Squared ' s2 10.2 'S ' s 10.2 'Lower CI ' lowerci 10.2 'Upper CI ' upperci 10.2;
title "Scheffe contrasts for this example. Refer to page 757";
run;
proc glm data=kenton;class design;
model cases=design/clparm alpha=.05;
means design/lines regwq;run;
 data Tire;
      input tire $ wear1 wear2 wear3 wear4 CellWgt;
      datalines;
tire1 24 25 25 25 1
tire2 27 25 25 25 1
tire3 25 26 25 26.5 1
tire4 26 23 28 23.5 1
;
   run;
/*pp. 760-761 give the ways to estimate the width of the Scheffe CI for planned comparisons.
  This is straightforward. I find it useful, though, to show power calculations for planned contrasts
  in GLM Power, as this is likely to be slightly more useful and you can adjust the experimentwise alpha
  in this if you wish. Four different scenarios are taken from Chapter 16, p. 721*/
proc glmpower;
   class tire;
      model wear1 wear2 wear3 wear4=tire;
      weight CellWgt;
      contrast "High vs. Low price" tire 1 -1 -1 1;
      contrast "Local vs. National" tire 1 1 -3 1;
	  contrast "All Pairwise" tire 1 -1 0 0,
	                          tire      1 0 -1 0,
							tire 	   1 0 0 -1,
							tire 	   0 1 -1 0,
							tire 	   0 1 0 -1,
							tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
         stddev = 2
         alpha  = 0.025
         ntotal = .
         power  = 0.8;
   run;
proc glmpower;
   class tire;
 model wear4=tire;
 weight CellWgt;
 contrast "High vs. Low price" tire 1 -1 -1 1;
 contrast "Local vs. National" tire 1 1 -3 1;
 contrast "All Pairwise" tire 1 -1 0 0,
                         tire      1 0 -1 0,
						 tire 	   1 0 0 -1,
						 tire 	   0 1 -1 0,
						 tire 	   0 1 0 -1,
						 tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
         stddev = 2
         alpha  = 0.025
         ntotal = .
         power  = 0.8;
 plot x=power min=.5 max=.95;
 title "power plot for Scenario 4";
   run;
proc glmpower;
 class tire;
 model wear4=tire;
 weight CellWgt;
 contrast "High vs. Low price" tire 1 -1 -1 1;
 contrast "Local vs. National" tire 1 1 -3 1;
 contrast "All Pairwise" tire 1 -1 0 0,
                         tire      1 0 -1 0,
						 tire 	   1 0 0 -1,
						 tire 	   0 1 -1 0,
						 tire 	   0 1 0 -1,
						 tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
 stddev = 2
 alpha  = 0.025
 ntotal = 40
 power  = .;
plot x=n min=24 max=480;
title "Sample Size plot for Scenario 4";
   run;
/*Assuming (p. 761) that the sample size for brand 4 is twice as large as the others*/
 data Tire;
      input tire $ wear1 wear2 wear3 wear4 CellWgt;
      datalines;
tire1 24 25 25 25 1
tire2 27 25 25 25 1
tire3 25 26 25 26.5 1
tire4 26 23 28 23.5 2
;
proc glmpower;
   class tire;
      model wear1 wear2 wear3 wear4=tire;
      weight CellWgt;
      contrast "High vs. Low price" tire 1 -1 -1 1;
      contrast "Local vs. National" tire 1 1 -3 1;
	  contrast "All Pairwise" tire 1 -1 0 0,
	                          tire      1 0 -1 0,
							tire 	   1 0 0 -1,
							tire 	   0 1 -1 0,
							tire 	   0 1 0 -1,
							tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
         stddev = 2
         alpha  = 0.025
         ntotal = .
         power  = 0.8;
   run;
proc glmpower;
   class tire;
 model wear4=tire;
 weight CellWgt;
 contrast "High vs. Low price" tire 1 -1 -1 1;
 contrast "Local vs. National" tire 1 1 -3 1;
 contrast "All Pairwise" tire 1 -1 0 0,
                         tire      1 0 -1 0,
						 tire 	   1 0 0 -1,
						 tire 	   0 1 -1 0,
						 tire 	   0 1 0 -1,
						 tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
         stddev = 2
         alpha  = 0.025
         ntotal = .
         power  = 0.8;
 plot x=power min=.5 max=.95;
 title "power plot for Scenario 4";
   run;
proc glmpower;
 class tire;
 model wear4=tire;
 weight CellWgt;
 contrast "High vs. Low price" tire 1 -1 -1 1;
 contrast "Local vs. National" tire 1 1 -3 1;
 contrast "All Pairwise" tire 1 -1 0 0,
                         tire      1 0 -1 0,
						 tire 	   1 0 0 -1,
						 tire 	   0 1 -1 0,
						 tire 	   0 1 0 -1,
						 tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
 stddev = 2
 alpha  = 0.025
 ntotal = 40
 power  = .;
plot x=n min=24 max=480;
title "Sample Size plot for Scenario 4";
   run;
PROC IMPORT OUT= WORK.Trainee 
            DATAFILE= "&location\Chapter17.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH17TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
data trainee;set trainee;
if training=1 then hours=6;else
if training=2 then hours=8;else
if training=3 then hours=10;else
if training=4 then hours=12;
proc print;
title "Piecework Trainees Example. Refer to Table 17.4, page 762";run;
proc sort data=trainee out=trainee;by training employee;
proc transpose data=trainee out=traineet prefix=employ;by training;
data traineet;set traineet;if _name_='product';
proc print data=traineet;
title "Information as it runs in Table 17.4, page 762, just for fun";run;
proc glm data=trainee;class hours;
model product=hours;
means hours/tukey alpha=.05 ;
title "Means and Anova summary table. Refer to Figure 17.6, page 763";run;
data trainee;set trainee;
chours=hours;
proc standard data=trainee out=trainee mean=0;var chours;
data trainee;set trainee;
chours2=chours*chours;
symbol1 i=rq v=dot c=black;
axis1 offset=(5);
axis2 order=35 to 65 by 5 offset=(5);
proc gplot data=trainee;
plot product*hours/haxis=axis1 vaxis=axis2 regeqn;
title "Scatterplot and fitted quadratic equation. Refer to Figure 17.7, page 764";run;
proc glm ;
model product=chours chours*chours;
title "Centered Quadratic Fit. Refer to Table 17.6a, Equation 17.51";run;
proc glm;class hours;
model product=hours;
title "Anova model. Refer to Table 17.6b, Equation 17.50";run;
proc glm;class hours;
model product=hours;
contrast "test of linear fit" hours -1.5 -.5 .5 1.5,
                              hours 1 -1 -1 1;run;
proc glm;class hours;
model product=chours chours2 hours;run;
title "Anova lack of fit test. Refer to Table 17.6c, page 765 and page 766";run;