%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter19;
PROC IMPORT OUT= WORK.castle 
            DATAFILE= "&location\chapter19.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH19TA07$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value widthfmt
 1='1=regular'
 2='2=wide';
value locfmt
 1='1=bottom'
 2='2=middle'
 3='3=top';
data castle;set castle;
format width widthfmt. location locfmt.;
proc means;class location width;var sales;
title 'Castle Bakery cell means Refer to Table 19.7 page 833';run;
proc means;class location ;var sales;
title 'main effect location marginal means Refer to Table 19.7 page 833';run;

proc sort data=castle out=castle;by width location;

proc boxplot;by width;
plot sales*location;
title 'boxplots of cells';run;
data castle;set castle;
locationf=location+.1*width;
proc sort data=castle out=castle;by locationf;
proc boxplot;
plot sales*locationf;
title 'boxplots of all cells';run;
proc glm;class location;
sales=location;
means location/hovtest=bf;
title 'one-way homogeniety of variance test- note in-class cautions about its use';
proc glm;class location width;
model sales=location|width/solution clparm;
means location|width;
lsmeans location/pdiff adjust=tukey cl ;
output out=castle r=resid p=predict rstudent=studdel h=leverage;
title "Castle Data Two Factor model. Refer to p. 837-840, p. 845";
title2 "Tables 19.9 and Figure 19.9, p. 842";
title3 "Table 19.10, p. 854";run;
symbol1 v=dot i=none;
axis1 offset=(10);
proc gplot;
plot resid*predict/noframe haxis=axis1;
title "Residual Plot. Refer to Figure 19.10, p. 843";run;
proc capability data=castle;var resid;
qqplot resid/noframe;
title 'probability plot using proc capability Refer Figure 19.10b, p. 841';run;
proc univariate;var resid leverage;
probplot resid/square normal noframe;
title 'probability plot using univariate';run;
axis1 offset=(10);
symbol1 v=none i=boxjt c=blue;
symbol2 v=none i=boxjt c=red;
proc gplot data=castle;
plot sales*location=width/overlay haxis=axis1 noframe;
title 'Line chart of mean differences by width. Refer to Figure 19.8, page 836';run;
data castle2;set castle;
if location=1 then d1=1;else d1=0;
if location=2 then d2=1;else d2=0;
width=2-width;
*width=2*(width-1.5);
d1width=d1*width;
d2width=d2*width;
proc means;var d1 d2 width d1width d2width;
title 'coding with dummy coding- replicates solution in glm';run;
proc reg;
model sales=d1 d2 width d1width d2width;
*loc1_2same: test d1=d2;
interaction: test d1width=d2width=0;
location: test d1=d2=0;
run;
data castle3;set castle;
if location=1 then d1=1;else d1=-2;
if location=2 then d2=1;else d2=-2;
*width=2-width;
width=2*(width-1.5);
d1width=d1*width;
d2width=d2*width;
proc means;var d1 d2 width d1width d2width;
title 'coding with effect coding- replicates statistical significance in glm';run;
proc reg;
model sales=d1 d2 width d1width d2width;
*loc1_2same: test d1=d2;
interaction: test d1width=d2width=0;
location: test d1=d2=0;
run;
proc gchart data=castle;
vbar location/sumvar=sales type=mean discrete noframe;
title "Vertical Bar Chart of means by Height. Refer to Figure 19.12, p. 851";run;
proc gchart data=castle;
vbar width/sumvar=sales type=mean discrete noframe;
title "Vertical Bar Chart of means by Width. Refer to Figure 19.12, p. 851";run;
proc glm;class location width;
model sales=location|width/solution clparm;
means location|width;
lsmeans location*width/pdiff adjust=bon alpha=.05 cl slice=width;
title "Examination of treatment means within level. Refer to p. 855";
title2 "This use of Bonferroni fudges the fact that you're interested in alpha=.10 for half of the ls means of interest";run;

/*Mathematics learning example. Nasty Hobbitses don't give us data, so we have to make our own*/
data learning;
do teaching=1 to 2;
do ability=1 to 3;
do i=1 to 21;
score=rannor(0);
output;
end;
end;
end;
proc standard mean=0 std=5.163977795 data=learning out=learning;var score;
proc means;var score;
proc format;
value abilfmt
1='1=Excellent'
2='2=Good'
3='3=Moderate';
value teachfmt
1='1=Abstract'
2='2=Standard';
data learning;set learning;
if teaching=1 and ability=1 then score=score+92;
if teaching=2 and ability=1 then score=score+90;
if teaching=1 and ability=2 then score=score+81;
if teaching=2 and ability=2 then score=score+86;
if teaching=1 and ability=3 then score=score+73;
if teaching=2 and ability=3 then score=score+82;
format ability abilfmt. teaching teachfmt.;
proc glm;class teaching ability;
model score=teaching|ability;
means teaching*ability;
title "Mathematics learning example. Refer (sort of) to table 19.11 p. 858";run;
symbol1 v=none i=boxjt c=blue;
symbol2 v=none i=boxjt c=red;
symbol3 v=none i=boxjt c=green;
axis1 offset=(20,20);
proc gplot;
plot score*ability=teaching/haxis=axis1 noframe;
title "Plot of estimated treatment means by method. Refer to Figure 19.13a, p. 859";run;
proc gplot;
plot score*teaching=ability/haxis=axis1 noframe;
title "Plot of estimated treatment means by abiilty. Refer to Figure 19.13b, p. 859";run;
proc glm;class teaching ability;
model score=teaching ability teaching*ability/clparm solution;
lsmeans teaching*ability/adjust=bon alpha=.975 cl pdiff ;
title "Bonferroni comparisons of differences. Refer to page 860";run;
proc glm;class teaching ability;
model score=teaching ability teaching*ability/clparm solution alpha=.05;
estimate 'mod greater than good gain?' teaching*ability 0 .5 -.5 0 -.5 .5;
title "Bonferroni comparisons of differences. Refer to page 860";run;
proc glm data=castle;class location width;
model sales=location width/solution clparm;
lsmeans location width/cl;
title "Pooled estimates. Refer to page 861";run;